text:
Processing proteomics data
[this question was also asked on BioStars]
In relevance to my research, I've been looking for proteomics data (control vs. diabetic) and I found a dataset in the article "Diabetes causes marked inhibition of mitochondrial
metabolism in pancreatic β-cells".
I am actually looking for the data that was used for creating figure 1. I checked the PRIDE dataset shared in the paper (https://www.ebi.ac.uk/pride/archive/projects/PXD012979) and I could find the raw files (RAW file, proprietary format) but I was not sure where to look for the processed data (does the .mzid file contain precessed data) and how to generate figure 1 (added below).

I checked a few workflows for processing the raw data (e.g. Progenesis and Knime- OpenMS). But the experimental conditions are missing in the PRIDE link.
I am not sure how to get started with this. Could someone please give some suggestions?
EDIT:
I understand from the answer posted below, one could use graph capture to digitize the values of relative abundance and create the heatmap. I am actually interested mainly in the enzymes involved in the glycolytic pathway Fig. 1 c,d of the paper (image shown above). Unfortunately, I couldn't find the relative abundance reported for proteins involved in this pathway.

You asked two questions

Is it possible to find from the dataset (i.e does the raw file contain information on the sample whether it is a diabetic sample or a control sample?).
Should the diabetic and control dataset be normalized separately?

Before I continue further I'd just like to point out the graph in the Supplementary information of the Nature research article of the paper in question. This graphic has not been edited in anyway.

I will not comment further but it would be very unusual for an analyst to produce a y-axis as done here. Just to labour/labor the point this is the appendix of a Nature research article (not 'even' a Nature letter).
The answer to the first question is you can't and will need to contact the corresponding author. There is a possible work around. The reason you need to contact the authors is:

The authors have not labeled the x-axis in any heatmap figure. This is not standard procedure and in this analysis you cannot link the heatmap values for that given mouse to the data base - because the authors have not declared which mouse is which. In addition, one cannot link the RNA-seq analysis to the proteomics data and they are paired samples, because only 3 diabetic mice are presented for the RNA-seq data but 4 mice for the proteomics data, therefore again it is not possible to link the two data sets (I understand this last point is not the precise question, but is important in the 'conclusion' below).
The link you provided had 12 samples - not 8 as in the paper.  It is very important to understand why there is a different number of samples (where one 'sample' is one mouse) in the proteomics data base compared with the apparent lower number samples presented in the paper.

To reiterate, only the authorship knows this information therefore they need to be contacted. The corresponding author is different from the data base curator who has been previously contacted (with no response). The corresponding author has to uphold the terms of publication as described at the end of the paper.

The workaround is to 'reverse engineer' the heatmap to its original values and is described in this article https://www.r-bloggers.com/2020/01/how-to-reverse-engineer-a-heat-map-into-its-underlying-values/ . This is only for R.

However, you still cannot link these values to the publicly available data sets, therefore until x-axis sampling is declared it is not a singular solution.
Question 2
Normalisation, I would like to point out in the appendix the authors transformed the RNA-seq data (log2) but did not do this for the proteomics data. For example aquaporins were upregulated 94 fold for proteomics data, but 4.61 fold using log2 RNA-seq data. The reason they did this is not clear and understanding this is quite important to the analysis. Why present an untransformed ratio in one analysis, but for the same experiment transform the values for precisely paired data?
The Appendix describes this information in the Supplementary Information in the following table.

Conclusion
This is a Nature paper and the authors are therefore under obligation, as per terms of publishing, to firstly supply all information on request and any information present in the published article must be explained in the Supplement Information. Those terms are stated at the end of the paper. It is clear the absence of labelling the x-axis in the heatmap in anyway (even if it is done in the legend), or the inability to link pairwise RNA-seq to proteomics data (which is very important), or the discrepancy in the samples on the database verses the samples in the paper, is not explained in Supplementary Information (Appendix).
Under these technical circumstances the authors are obliged to supply the missing information on request. You may require a chain of request escalation within your organisation because they may not respond to the initial requests, but you don't know unless you try. However, at some level in the organisation, e.g. Director, the terms of publication would be upheld by authorship in question.

First answer

But the experimental conditions are missing in the PRIDE link

The experimental conditions of the heatmap are as follows: each box is a separate animal and the heat map colors (colours) are a log2 representation of Figure A (below). In other words to understand the heat map using Figure A is the starting point. You can obtain the values from this analysis easily using a graph capture package using Figure A (not Figure B). Its a long while since I've used graph capture, so I forget the packages. If you feel there are anomalies you can then contact the authors for the precise values.
log2 is a bit weird, natural log would have been preferred:
FIGURE A
"Abundance of the indicated proteins involved in lipid synthesis, measured by proteomics, in islets isolated from control (black, Ctrl, n=4) and 2-week diabetic βV59M (white, Diab, n=4) mice. Each data point indicates a separate mouse. Mean±S.E.M. **p<0.01, ***p<0.001. HMGCS1, 3-Hydroxy-3-methylglutaryl-CoA synthase 1; HMGCR, Hydroxy-3-methylglutaryl- CoA reductase; DHCR7, 7-dehydrocholesterol reductase; ACLY, ATP citrate lyase; GPD1, glycerol phosphate dehydrogenase; ACC1, acetyl-CoA carboxylase 1; AACS, acetoacetyl- CoA synthetase; FASN, Fatty acid synthase."
FIGURE B
"Heat maps of mRNA and protein expression of the indicated lipid metabolism genes in islets isolated from control and 2-week diabetic bV59M mice. Each box corresponds to a different animal. Colour indicates log2 fold-change."
https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-019-10189-x/MediaObjects/41467_2019_10189_MOESM1_ESM.pdf

Data availability
RNA data
https://www.ebi.ac.uk/ena/browser/view/PRJEB31793?show=reads
Checks out ok.
Proteomics
I checked one link only ...
http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD012979
I agree something doesn't look right here just in the way it is reported, because the paper should highly polished and they did miss the link out. This could be word limit issues (Data availability might be included). The curator is listed as roman.fischer@ndm.ox.ac.uk and that would be the person to contact.

Thanks so much. Please check my edit. I actually wrote to the curator but couldn't get clarification. I was suggested to process the raw data. I'm not sure if you had a chance to look at this https://www.ebi.ac.uk/pride/archive/projects/PXD012979. This has the raw data but I am not sure out of the 8 samples which corresponds to control/diabetic. Could you please have a look at this link?

Firstly, please give the curator time to respond, at present they will be inundated and it will be one person supporting loads of stuff. Anyway many thanks for your update I now understand, thanks. I will have a look in more depth this week @Natasha. My gut feeling is, if it was easy you would have figured this out.

"The link you provided had 12 samples - not 8 as in the paper." I found 12 + 8 samples on PRIDE (12 files with file name prefixed with
FL0461_MSQ969b_FAshcroft and 8 files prefixed with FL0408_MSQ911_MariaRohm). But I am not sure whether these 8 samples prefixed with FL0408_MSQ911_MariaRohm correspond to the 8 samples reported in the paper. Could you please comment?

I'm afraid there are definitely 12 samples. The "8" are not samples but files, two files are identical one is zipped one isn't, the "...all.mgf" and "...all.mgf.zip" and two corresponding to the same label "..._M1.raw". The 12 samples  are labeled "...sampleA.RAW", "...sampleB.RAW" ... to "...sampleK.RAW" "...sampleL.RAW". The intervening files are listed in alphabetical order.

In summary, unfortunately the 12 samples in question are files singularly labeled "sample" from A to L. I would advise clearing this up with the authorship, they are the only people who know.

I wrote more than 10 emails to 4 authors of the paper. But no one clarified about the samples. I was told to do a PCA and find out which sample corresponds to what :/. Reg. files, "FL0408_MSQ911_MariaRohm_544_F31.raw", "FL0408_MSQ911_MariaRohm_544_M1.raw" etc also contain the raw data.  ( "...all.mgf" and "...all.mgf.zip"  were identical, I agree).

This would be expected. Thus I'm suggesting you use escalation, i.e. for the next person above you to inquire to the corresponding author. What I am saying is the authors have to deliver the information because that is the terms of publication in Nature. There will be a tipping point where the seniority of the individual will force the authors to comply. If they did not the journal will enforce its publication policy, which ultimately would lead to an enforced official 'correction'. The authors will want to avoid this scenario at all costs.

The response above contains the details regarding the technical points that define the authors responsibilities of course in which your inquiry is part. If you simply say 'the database is jumbled' that would not really work. If you are still stuck you simply show this response to the person immediately above you, i.e. your immediate senior, together with the formal terms of publication (in the journal). They should understand how to proceed.

could you please have a look at this post https://bioinformatics.stackexchange.com/questions/18928/heat-map-of-protein-expression-from-normalized-abundance

I got a reply yesterday from the corresponding author on the diabetic vs. control labels of the samples. Just wanted to post you on this.

Thanks for the update. They will likely have seen this post by now and be aware of the discussions therein.

I'm glad I opened a discussion here. Thanks a lot for your support. I tried processing the raw files (8 RAW files with labels MariaRohm) using MaxQuant and got the following result (file). But the magnitude of LFQ intensities (I hope this is the normalized abundance, please correct me if I am wrong) differ from the values shown in the paper. I tried to compare it for PDK1 panel e in Fig2 of the paper.

I updated this above comment in a new post (https://bioinformatics.stackexchange.com/questions/18929/setting-up-maxquant-parameters-for-processing-proteomics-data). Thank you

Discussions should be carried out in chat. Upvoting comments means that discussions get confusing, because comments will appear out of order. Where possible, please edit your questions and answers, then delete the discussion comments.


text:
Multiple file Python script
I am currently doing an analysis requiring a multiple-input script. I have a directory filled with 900+ samples.
Each sample is comprised of two-files in a fastq format. This is an example of a sample:
MIP_REV_BC_01_01A_S1_R1_001.fastq.gz
MIP_REV_BC_01_01A_S1_R2_001.fastq.gz

The script has a Read1 & Read2 variable, and I need to make it so the R1 file is inputted for Read1 and the R2 file is inputted for Read2 for every sample. To be clear, the script is supposed to be executed on every sample separately.
This was an attempt:
with open(filename) as file:
    lines=file.readlines()
    for line in lines:
        parts[]=line.split("_")
        if parts[7]="R1":
            r1[].append(line)
        else if parts[7]="R2":
             r2[].append(line)

We were bogged down by syntax errors, but either way, this script doesn't seem to do what we wanted it to. However, it captures the idea of what we're trying to achieve.
Question I want to input a directory into the function so it automatically runs the script on the samples in the directory.
Any idea what I should do?

is the idea that you want to process both files at the same time and not store them in memory?

Edited with explanation below. I don't really understand the objective, because I initially thought you were editing the content of fastaq file, but you seem to be editing a file which contains fastaq filenames. Is that correct? I know the syntax but I don't know really what the goal is.
Note I'm using the pathlib which has a glob function, do consider rglob under pathlib.
import re
import pathlib as Path

# pathlib is super, good to habit to use it
mypath = Path('/Users/fewbugs/workingd')
allfilenames = mypath.glob("*.fastaq.gz")
r1, r2 = [], []
x = re.compile('_')
for filename in allfilenames:
    # the 'r' is super important
    with open(filename, 'r') as file:
    # there is no need for a readlines()
        for line in file:
    # split via re, it is better to compile (faster)
        parts=x.split(line)
    # you must use '=='. REMEMBER lists start at 0 not 1
           if parts[6]=="R1":
    # if you want to append to a list, declare the list first
            r1.append(line)
    # there is no such thing as 'else if' REMEMBER lists start at 0 not 1
            elif parts[6]=="R2":
                 r2.append(line)
    # when testing good idea to trap the stuff you're missing
            else:
                 print ("This one escaped %s" % (line))

Please be CAREFUL with my spacing, everything needs to be multiples of 4. I've not checked the code.

Forget this bit I thought you were processing the content of the fastaq file.
outfile1 = Path(mypath.parent, 'outfile1')
outfile2 = Path(mypath.parent, 'outfile2')

with open (outfile1, 'w') as f:
    for item1 in r1:
        f.write("%s\n" % item1)
with open (outfile2, 'w') as f2:
    for item2 in r2:
        f2.write("%s\n" % item2)

A simpler way?
... it looks like you're just processing filenames...
Might I suggest just globing no need for all this stuff.
mypath = Path('/Users/fewbugs/workingd')
r1filenames = mypath.glob("*_R1_*.fastaq.gz")
r2filenames = mypath.glob("*_R2_*.fastaq.gz")

If that is the aim, thats it! Just be careful about the glob versus rglob.
Idea If you want two paired lists this might be better done in a nested comprehension, e.g. a counter could be inserted at the central glob (*).
That may very well be what you're after. i.e.
R1 = ['stuff_R1_001.fastaq.gz', 'stuff_R1_002.fastaq.gz' 'stuff_R1_003.fastaq.gz']
R2 = ['stuff_R2_001.fastaq.gz', 'stuff_R2_002.fastaq.gz' 'stuff_R2_003.fastaq.gz']

source: Stack Exchange

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

sort: votes


text:
How can I separate upregulated genes from downregulated on a boxplot?
I have a table of fold change shifts for two different groups. Some of them are upregulated, others are downregulated. I wonder if it is ok to subset only positive values and to boxplot them to show only upregulated genes and vice versa?
     Gene   REDmean   GRNmean
1    gen1   4.918239  4.547266
2    gen2   3.635476  3.447484
3    gen3   -3.759263 -3.283889
4    gen4   3.962326  2.818917
5    gen5   -3.474217 -3.123105
6    gen6   3.653160  2.907355

Try again, use separate boxplot columns i.e. upregulated vs. downregulatated that are created using a dataframe e.g. pandas. For pandas use a filter to create a new column of 'UpDown', i.e. a boolean of upregulated or downregulated, and then use the groupby('UpDown') command. Send that into matplotlib boxplot function and it should work. The 'Gene' column needs removing (I think). My pandas is too rusty to knock out the code. Someone here might show a solution in R, but if you attempt a pandas approach and post it, someone will probably run with it.

Basically you'll reorder and pivot the data using a dataframe using the up- down-regulated filter (I'm not sure if 'pivot' is a correct term). What I don't like is your use of the word 'mean', for the real thing you need the raw data ... but this would be cool for creating the filter. "'Gene' column needs removing" well thats a bad idea (thats your results), but for the plot it will cause it to fall over if its dragged in on a dataframe - thats what I meant in the comment above.

Thank you. As I have asked in my previous posts I am still thinking if this is a valid approach to compare the mean fold change between two groups. And once again, I wonder if I can compare their fold changes (FC), say, for gene1 in RED group I have 5 for control and 25 for treated sample and in GRN group (for the same gene1) I do have 2 for control and 10 for treated. For REDs I do have a 5X FC increase (from 5 to 25) and for GRN it is also 5X (from 2 to 10). For me there is no difference between RED and GRN groups based on gene1.

The idea of the 'REDmean' and 'GRNmean' is that I have multiple samples, so I have taken a mean control vs treatment ratio (that are presented in a table above). The question is if I have the right to do that, and if yes - can I distinguish them just by separating negative FC from positive?

AS A PLOT - I have no problem. As a statistical test ... I think we covered the stats and I personally would not favour change-fold ... but if is an approach to standardisation (its a separate question) that is different. However, others might agree with you, i.e. a test on change-fold. The operation you are describing here to manipulate and plot is a serious question that requires formal code in R or Python. "Try again" ... sorry this was a personal comment, my first comment got deleted so I "tried again", this could be misinterpreted however.

It is good to keep questions distinct, because this question over a very large dataset isn't trivial, i.e. to provide the code that will perform the pivot. I honestly can't remember fully how to do it without looking it up (I've done similar things in the past).

I'll be working on pandas again shortly (couple of weeks). I'll be able to answer the question then if someone doesn't do this before.

Great! Thanks a lot for the help

Let us continue this discussion in chat.

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


text:
Locate the saved template of minknow
I am using minknow on ubuntu. After saving the template (settings), where is that file locally stored in the system?
Also what is its format/extension?

It would help to understand a little better what you are actually doing. I can find a github page for MinKNOW and a page with a ton of tutorials. Where have you looked to try to find this information? Have you consulted the reference manual?

Hi! Thanks for responding. I want to share those files across a few systems, automate and load the exact same settings.
There is no content regarding the local storage for the templates.


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

There is no obstacle in performing a DE analysis and all downstream on a 6 vs 2 comparison. The only thing to keep in mind is that the representation of the control is not particularily good, but you can technically perform the exact same analysis as if it was 6 vs 6, e.g. with DESeq2, edgeR or limma-voom/trend.

Agree. In the absence of evidence for variation within a single sample type, DESeq2 can make guesses at variation based on overall variation for all samples.


text:
Seurat Integration
I'm following the instructions for integration:
https://satijalab.org/seurat/articles/integration_introduction.html
And it's taking a while to run:
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
    anchor.features = features)

I'm integrating 6 datasets.
Is this function designed well for large amount of data?

Are you using the same data from the integration tutorial or using your own data? How many cells/genes are in your dataset?
Seurat's default integration method (CCA) is known to be runtime/memory intensive. It can handle large datasets but may require lots of CPU cores/memory. I have found that the amount of RAM also seems to increase as more CPU cores are used. That said, I have been able to run the integration tutorial with the included dataset on an 8GB laptop.
The easiest solution is to simply use more resources when integrating: increase the amount of RAM and consider using multiple cores (there is a separate Seurat tutorial on parallelization). If you are encountering runtime/memory issues and you cannot scale up your resources, I would suggest looking at alternative integration methods:

Seurat includes another integration approach, RPCA which is supposed to be less memory intensive at the cost of being more conservative with integration
There are several other integration methods that provide similar results. I have used Harmony and found that it requires significantly less resources than the Seurat integration approaches (both CCA and RPCA). A quick google search of single cell RNA-seq integration methods will turn up other popular methods with benchmarks to compare and contrasts different approaches.


text:
How will Seurat handle pre-normalized and pre-scaled data?
I want to use data from the following dataset:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84371
The data has been TPM normalized, which is not ideal for clustering but I have to work with what I have. TPM is a bad normalization method and it should not be used for these analyses because its laden with a lot of assumptions.
Presumably it has already been scaled, can I simply skip these steps and then run FindVariableFeatures() followed by the rest of the pipeline from the clustering tutorial or do the methods used assume Seurat specific preprocessing and using TPM will screw them up?
Unfortunately in this case TPM is not reversable.  https://btep.ccr.cancer.gov/question/faq/what-is-the-difference-between-rpkm-fpkm-and-tpm/#:~:text=Here's%20how%20you%20calculate%20TPM,divide%20this%20number%20by%201%2C000%2C000.
TPM is calculated as a function of the kilobases of the gene in the reference genome they used and the reads per gene across all samples (aka all cells).  They have not specified which reference genome was used in the main article, and they do not supply the read counts in the metadata or the main file.  So in this case the TPM is irreversible.
I think I can get raw reads from SRA and do alignment myself. It's annoying that raw reads were not provided but as someone who's lab entirely does computational analysis on open datasets I can tell you that this is not the most egregious thing I've ever seen. As I said before it would be easier for me to get the fastq files and do my own alignment which is still a ton of work the main authors could have prevented by simply posting the read counts.
Clustering Tutorial:
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

TPM means transcripts per million. It is obtained by dividing the read counts by the length of each gene in kilobases, resulting in reads per kilobase (RPK). All the RPK values in a sample are counted and then divided by 1,000,000.

I don't do transcriptome analysis, it ain't my thing, however I do understand statistical analysis as well as the underlying issue regarding the public availability of molecular data ... I agree with the OP its not ideal. However, yes the OP can continue with 'clustering', personally I definitely prefer it if the transformation was removed, or analysed from the raw data.
Summary The authors of the data set are obliged to provide the raw data, presuming its published. Why? This is because whilst the OP can perform the analysis - it is very likely the OP will be asked to repeat the analysis on the untransformed data (again note this ain't my area, so 'very likely' might be 'possibly' in this field). At this point the original authors have to declare their untransformed data, if this was a formal request to the OP's analysis.
Please read the caveats and discussion.

Firstly could the OP reverse the normalisation? Its just a transformation applied to all the data which in their analysis resulted in a normal distribution, I don't know TPM but if this was log-transform you would supply 'un-log' it, square-root then square it, reciprocal do another reciprocal will 'de-transformation'. Can TPM be undone? I agree it ain't ideal, from several perspectives but it gets things moving quickly.

Ask the authors for the raw data - this is what I would do. They are obliged to provide it if the data is published. Virtually all journals have this policy in place, however I agree that ASM remain an outlier in this respect (who are significant publisher). Any other journal will return to the authors and point out their 'check-boxes' when they submitted the manuscript. It is etiquette to go to the authors first and provide sufficient chance for them to respond, prior step 3 (below) - giving them the absolute maximum benefit of the doubt.

What you do if the authors don't respond? The OP would then contact the journal explaining: a) why you need the data, i.e. the analysis performed; b) explain the analytical significance of the raw data, i.e. arguments the OP has raise 'not ideal' and general concerns raised here. I do have experience in this regard and the authors will oblige if asked by the journal.

Searaut's cluster analysis certainly used tSNE ... I believe it uses dimensionality reduction such principle component analysis (PCA) - I don't know which they use - and then tSNE analysis (which I think you are calling 'clustering*'). This is a very robust 'clustering'* (please see footnote), there is no question this is as good as it gets. If that is singly your goal I agree that dimensionality reduction calculations don't need transformations, transformations are needed for parametric statistics and other areas of model fitting.

The problems with transformation
I agree the authors should not have done this because:

Transformations are not perfect, if its wrong - you're stuck.
Not only do you not need to transform the data, its better you don't - if untransformed data results in nice 'clusters' thats the preferred result. Dimensionality reduction + tSNE gives great results thats well known.

Corroborating the OP's concern
Dimensionality reduction calculations do not need a normal distribution because its sectioning up the observable variance to maximise the classification of the data. My opinion - in statical analysis - is if you need to transform, THEN perform dimensionality reduction (PCA-stuff) THEN tSNE, that data has been heavily manipulated. It doesn't mean to say the conclusions are not valid but its been put through a lot and would STRONGLY SUGGEST THE RAW DATA HAS NOT GIVEN GOOD RESULTS. Thus, personally I would be cautious, why should be OP be put under unnecessary caution?
It is important to note this is not my field and the field will have its own standards. If transformation is routine for all statical calculations (except non-parametric where it may not be valid) then authors of the data will argue they are correct and no further data needs to be released. Personally, I really doubt it, but we'll see what others say.
Footnote * clustering is a statistical method usually resulting in a tree, it is generally avoided because there are so many different clustering methods giving different results and no agreement on which is the 'optimal' method. I think when the OP say 'clustering' they refer to clusters following a tSNE analysis, which will be a bivariate plot.

From the comments Failure to declare a reference genome is a formal omission that has to be rectified and the journal is obliged to do this. Thus the OP would simply contact the journal, they may tell the OP to contact the authors and when the authors respond the response is passed back to the journal. If its an online journal you can simply leave a comment in the 'comments' section: the authors will answer. In this instance officiating via the journal is preferred to ensure the information is verified by all parties - I presume the OP does not know the authors. From the comments I get the gist the OP is working via supervision, and would therefore just notify  their supervisor who should officiate on OP's behalf (they might then delegate it back to the OP - but at least the supervisor is formally aware of this).
In context, stuff happens and there is formal process to ensure parity is achieved its simply patience in working through it.

Down vote In response to the downvote ... unpopularity happens in data set transparency its part of the process ... I was (at least in part) attributed to changing journal group policy on data availability, stated from within the journal. Maybe they joking - but they certainly all knew about what was going on. The specific issue was a prolonged dialogue over access to a specific data set. It does demonstrate that most journals are concerned if an author does not cooperate with data transparency, in fact the journal in question explicitly said this and they were not joking.

With regards to your comments about the dimensionality redcution this is par for the course in RNA-seq one point to note though is that the results of the clustering should be robust across different embeddings (tSNE, UMAP, spectral clustering, CCA, PCA etc).

Thanks for the update. Regarding the last point ... the authors of the data don't have to be cooperative, but they have to be transparent. What you have said is there is a chance they have not been transparent and that is rarely accepted with molecular data.

Seurat doesn't use tSNE for clustering; tSNE is used only for visualisation.

I notice that the paper points to the Broad website for data accession; have you looked there?
https://singlecell.broadinstitute.org/single_cell/study/SCP1/-single-nucleus-rna-seq-of-cell-diversity-in-the-adult-mouse-hippocampus-snuc-seq
For specifics on the normalisation method, you could look at their associated github repository:
https://github.com/yinqingl/nucseq_analysis
Given that the paper was published in 2016, they are TPM counts, and you have access to the raw reads, my recommendation would be to regenerate gene count tables from those raw reads. I expect it will be more time consuming, less accurate, and require more effort to track down the gene count table than to create one from scratch. You can use Salmon Alevin to do this fairly quickly:
https://salmon.readthedocs.io/en/latest/alevin.html
[bearing in mind that there may be some footwork required to work out how to get Salmon to interpret the barcoding method that was used]


text:
Is there a generally accepted method used to impute missing DNA methylation data (probe beta-values)?
I'm looking at DNA methylation (DNAm) data such as TCGA (e.g., BRCA, KIRP, KIRC, etc.). Currently trying to build use my model to predict DNAm age on test sets, but many of the data sets  are missing key probe values used in my clock/model.
I have inspected approximately 30 GEO data series and the following TCGA cancer samples: KIRP, KIRC, LUAD, LUSC, BRCA, THCA. They are each missing key probes I am using.
Is there a common way to impute the missing values in R?

Could you please clarify what "my clock/model" means, and what specifically is missing? How do you know that probes are missing? If this is from methylation chip data, have you confirmed that the probes exist in the probes used for that chip?

Hi, could you show what kind of data do you have and how you identify these missing probes? This might help having some code to build upon and answer your question. I think there isn't any general way to impute missing values, so knowing the data and its sparsity might help proposing imputation methods.


text:
Error occurrence after merging files with bcftools: wrong number of fields?
This question was also asked on Biostars
I have multiple vcf of CASES and CONTROLS variations annotated by VEP, SNPEff and SnpSift.

First pair vcf -> only variations| CASES and CONTROLS
Second pair vcf -> variations + SnpEff | CASES and CONTROLS
Third pair vcf-> variations + SnpEff+VEP+SnpSIFT| CASES and CONTROLS | at so on
The data has non-overlapping samples. I tried to merge corresponding vcf of CASE and CONTROLS into one file with bcftools merge.

The code is:
 bgzip -c file.vcf > file.vcf.gz
 tabix -p vcf file.vcf.gz
 bcftools merge -o merged.vcf.gz -Oz file.vcf.gz file_2.vcf.gz

Unfortunatelly in files with SnpSift annotations i got this error:

Error at chr1:20724637: wrong number of fields in dbNSFP_ada_score?

Please anybody knows how to solved this error or what does it mean?
grep header:
INFO=<ID=dbNSFP_ada_score,Number=A,Type=String,Description="Field 'ada_score' from dbNSFP">

-G on chr1:
CHROM POS ID REF ALT QUAL FILTER INFO
chr1 20724637 rs1280663306;rs1445195248 TGGGAGGGAGGGAGAGAGGT TGGGAGGGAGAGAGGT,TGGGAGGGAGGGAGGGAGGT,TGGGAGGGAGGGAGAGAGGG,TGGGGGGGAGGGAGAGAGGT,TGGGTGGGAGGGAGAGAGGT,TGGGAGAGAGGGAGAGAGGT 2.42927 . ANN=TGGGGGGGAGGGAGAGAGGT|splice_region_variant&intron_variant|LOW|SH2D5|SH2D5...SNV|1||||||||||||||||||||||||||||||||||||||||||||;dbNSFP_ada_score=1.43020658327008E-4,0.00157457299236331,5.16496450143413E-5;dbNSFP_rf_score=0.016,0.054,0.0;dbNSFP_PHRED=5.632,0.077,6.365,7.459,2.716
chr1 20724637 rs1280663306;rs1445195248 TGGGAGGGAGGGAGAGAGGT TGGGAGGGAGAGAGGT,TGGGAGGGAGGGAGGGAGGT,TGGGAGGGAGGGAGAGAGGG,TGGGGGGGAGGGAGAGAGGT,TGGGTGGGAGGGAGAGAGGT,TGGGAGAGAGGGAGAGAGGT 2.42927 . ANN=TGGGGGGGAGGGAGAGAGGT|splice_region_variant&intron_variant|LOW|SH2D5|SH2D5|transcript|NM_001103161.2|protein_coding|5/9|...||||||||||||||||||||||||||;dbNSFP_ada_score=1.43020658327008E-4,0.00157457299236331,5.16496450143413E-5;dbNSFP_rf_score=0.016,0.054,0.0;dbNSFP_PHRED=5.632,0.077,6.365,7.459,2.716

dbNSFP_ada_score is an ensemble prediction score, as calculated by dbNSFP.


text:
Number of collinear genes are zero while running MCScanX
I am trying to do a synteny analysis using MCScanX. As per requirement I have done all by all blast between two sets of protein from two closely related species. And also prepared the gff file. The file formats are as follows:
cyp_dan.blast
NP_001347641.1  XP_042591522.1  78.31   166     32      1       1       166     1       162     5e-91   267
NP_001347641.1  XP_042577882.1  60.24   166     62      2       1       166     1       162     4e-63   196
.
.
XP_704272.2     XP_018953536.1  31.03   87      59      1       275     360     90      176     8e-07   50.4
XP_704272.2     NP_001005958.2  33.33   63      38      1       273     331     38      100     1e-06   48.9
XP_704272.2     XP_018969808.1  31.33   83      56      1       275     356     91      173     2e-06   49.3```

cyp_dan.gff
NC_056572.1     gene-LOC1221465.1       2966    7093
NC_056572.1     id-LOC122146575.1       2966    3166
NC_056572.1     id-LOC122146575.1       3246    3357
.
.
NC_002333.2     gene-trnP;Dbxre.1       16527   16596
NC_002333.2     rna-trnP;Parent.1       16527   16596
NC_002333.2     exon-trnP-1;Par.1       16527   16596

cmdline: /home/software/MCScanX/MCScanX cyp_dan
Output:
############### Parameters ###############
# MATCH_SCORE: 50
# MATCH_SIZE: 5
# GAP_PENALTY: -1
# OVERLAP_WINDOW: 5
# E_VALUE: 1e-05
# MAX GAPS: 25
############### Statistics ###############
# Number of collinear genes: 0, Percentage: 0.00
# Number of all genes: 108556
##########################################

Here, both of the species are from same family since there is no chance to have "zero collinear gene".

I'm not familiar with the MCScanX tool, but I see two potential issues that are worth investigating.

The cyp_dan.gff file does not look like GFF or any of its common variants (GFF3, GTF, etc.). Again, I'm unfamiliar with MCScanX, so I don't know precisely how its wants the input data formatted. But if it wants a GFF file, I would double-check to make sure that what you have matches what the tools expects/requires.

The gene/protein IDs in the cyp_dan.gff file do not seem to match the query or subject IDs in the cyp_dan.blast file. No tool is going to be able to find collinear adjacent genes if the IDs in the BLAST results don't match the IDs in the gene annotations.

Hi,
Thanks for your response.

Original gff format is not allowed to use. It requires modification with having chr name, geneID, start and end position only.
I took out the necessary info using following cmd:
awk -F "\t" '{a=substr($9,4,15)}$0~/gene/{print $1 "\t" a".1" "\t" $4 "\t" $5}' A.gff > A_new.gff

I have extracted the info from the NCBI annotated file and didnt make any change in ID info by my self. And after extracting this is how it looks. Can you please give me any suggestion to find exact match between blast and gff file?


text:
De novo antibody sequencing fromMS/MS Ion Trap proteomics raw signal
I've read about different proteomics tools that can do de novo assembly of MS/MS Ion Trap proteomics signatures for de novo antibody sequencing. Given the way CDR3s of antibodies are created, together with Somatic Hypermutation, can generate a huge diversity in the CDR regions of the heavy and light chains, the reference-based methods can fall short of determining the full sequence of an unknown antibody without a set of reference sequences including those with its CDRs.
Are there any other tools that can determine antibody sequences from MS/MS proteomics raw signals?

A few notes on data generation first:
There are currently two broad approaches for the generation of bottom-up or "shotgun" MS proteomics: data-dependent acquisition (DDA) and data-independent acquisition (DIA).1 In tandem MS (MS/MS), the DDA approach only puts forward certain peptides generated during the first cycle of MS for fragmentation during the second cycle, while with the DIA approach, all peptides generated during the first MS cycle can be fragmented in the second round.
The recently-introduced Orbitrap Fusion mass spectrometry permits various types of MS2 acquisition methods. There are four main MS2 strategies:

HCD-OT (higher-energy-collisional-dissociation with Orbitrap detection),

HCD-IT (HCD with ion trap, IT),

CID-IT (collision-induced-dissociation with IT) and

CID-OT on Orbitrap Fusion

Reading from the biorxiv preprint of Stitch:
https://www.biorxiv.org/content/10.1101/2022.03.07.483237v1
"With the aid of specialized software packages like DiPS, Supernovo, or ALPS, full and accurate sequences of the heavy and light chains can be reconstructed with the MS/MS spectra of the digested peptides."
Before running Stitch or similar tools, a peptide profiling tool needs to be executed on the .raw or .mgf data alone.
STITCH: https://www.biorxiv.org/content/10.1101/2022.03.07.483237v2
This currently requires the output of analysing .raw data via PEAKS, which is not free unfortunately. Authors currently looking into using other free de novo sequencing software for input into Stitch. The software also allows for FASTA format, which can then be executed on the output of PepNovo, Novor or similar.
DiPS: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5461544/
ALPS: https://www.nature.com/articles/srep31730
"The results show that ALPS was able to assemble three complete monoclonal antibody sequences of length 216–441 AA, at 100% coverage and 96.64–100% accuracy."
DeepNovo: https://www.pnas.org/doi/10.1073/pnas.1705691114
"This field has been actively studied over the past 20 y, and many de novo sequencing tools have been proposed, such as PepNovo, PEAKS, NovoHMM, MSNovo, pNovo, UniNovo, and Novor among others (6–19). The recent "gold rush" into mAbs has undoubtedly elevated the application of de novo sequencing to a new horizon (20–23)"
"DeepNovo considerably outperformed state of the art methods, achieving 7.7–22.9% higher accuracy at the amino acid level and 38.1–64.0% higher accuracy at the peptide level"
Available for download, requires activation via sending email to author, which then points to the company behind PEAKS, which will give you a 15 day evaluation license.
NovoHMM: the download page leads to a missing URL
DeepResolution2: https://www.sciencedirect.com/science/article/abs/pii/S0039914022002119
Capable of analyzing GC-MS data rather than MS/MS Ion Trap
DePS: https://arxiv.org/abs/2203.08820
DeepNovoV2: https://github.com/volpato30/DeepNovoV2
DPST: https://arxiv.org/abs/2203.13132
Link does not work https://github.com/Yan98/DPST
PepNovo: http://proteomics.ucsd.edu/software-tools/531-2/
The new version is from 2012, includes models for CID and LTQ, but not HCD models.
pNovo: http://pfind.org/
Includes models for CID, HCD, ETD and HCD+ETD.
Novor: https://www.rapidnovor.com/
Filling in a form, someone will contact you. Doesn't seem to be free.
METHODOLOGY used so far:
1 - ThermoRawFileParser: Convert .raw data to .gcf
2 - PepNovo: produce list of peptides, one-liner to convert to .fasta
3 - Stitch: analyse given a list of V(D)J+C references from templates and the .fasta file from above.


text:
How can I add a ";" at the end of each fasta header line in a file?
I have a large fasta file with RNA sequences. I need to add a ; at the end of each header line (lines beginning with >) in the file.
For example, I have:
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

And, my desired output is:
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

One way, using sed and a regex address to select only the header lines and apply a substitution to append a semicolon to the end of each line:
sed '/^>/s/$/;/' in.fa > out.fa

Or, to edit the file in place:
sed -i '/^>/s/$/;/' in.fa

Or one way using awk and a ternary construct to print the line followed by a semicolon or an empty string, depending on whether or not the line is a header line:
awk '{ print $0 (/^>/ ? ";" : "") }' in.fa > out.fa

Alternatively, you may prefer:
awk '/^>/ { print $0 ";"; next } { print $0 }' in.fa > out.fa

Which is the same as:
awk '/^>/ { print $0 ";"; next }1' in.fa > out.fa

If you have perl installed, all you need is:
perl -pe 's/$/;/ if /^>/' in.fa > out.fa

... using the legendary Perl pie
perl -p -i -e 's/^(>.*)/$1;/' mybacteria.fa

Input
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG
A

Output
>Bacteria;Cyanobacteriales;Phormidiaceae;Lyngbya PCC-7419;Lyngbya sp. PCC 8106;
GGAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUCUGCUUAACACAUGCAAGUCGAACGGACUCUUCGGAGUUAGU
>Bacteria;Cyanobacteriales;Nostocaceae;Nodularia PCC-9350;Nodularia spumigena CCY9414;
GAGAGUUUGAUCCUGGCUCAGGAUGAACGCUGGCGGUAUGCUUAACACAUGCAAGUCGAACGGUCUCUUCGGAGAUAGUG

Details I copy pasted your first code into a file mybacteria.fa then ran the Perl pie script. Perl pie is very quick, there alternative regexs (regular expressions) that would be slightly more efficient, but Perl pie is fast regardless. In fact you don't need to use regex in this example and would be faster via a if $line == '>'; declaration (certainly true for Pythonic equivalents if line == '>') but it doesn't offer the ease and flexibility, which is combined with speed, that Perl pie does. The code will work in Linux (any flavour/flavor) or Mac OSX Darwin (i.e. Terminal app)

It will work on any flavour/flavor so long as Perl is installed. Note that not all distributions ship with Perl pre-installed. Also, the capture group is not necessary - all that's needed here is: perl -pe 's/$/;/ if /^>/' in.fa > out.fa. Note that using -i without an extension will modify the input file in-place without making a backup of the original.

In fairness I did point that out in the discussion below the code. The efficiencies are generally not noticeable. However, most regex manipulations take place within the string rather than at the end requested here, and capture can be one way to deal with different situations beyond this one.

Fair point. Perl is included in most popular Linux versions including Ubuntu, Debian and Suse. Perl isn't preinstalled on RedHat or freeBSD after release 5 and nobody uses Kali in bioinformatics. Let me know if there's anything missed.

Thank you very much -- this worked well and was simple to execute.

@CallahanMcGovern if an answer worked for you, can you please mark it as accepted?

I recommend to use Biopython as a clearer way.
from Bio import SeqIO

records = []
for record in SeqIO.parse("file.fas", "fasta"):
    record.description += ";"
    records.append(record)
SeqIO.write(records, "file.fas", "fasta")

Before:
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA
TTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGATCCAATAATCCTGCTGGAAT
TACATCAGAATCTGACAAAATTCCATTTCACCCATACTACACCATTAAAGACATTCTGGGAATCCTATTC
CTTATTTTTATACTTATAGGTTTAGTGCTATTTACGCCAGACCTACTAGGAGACCCAGATAACTACACAC
CCGCCAACCCCTTAAGCACCCCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCT
ACGTTCCATTCCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCATTCTGATTCC
TAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAATATCCCTTCATCATCATCGG
CCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTCCT

>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAGATCGCAACAG
GTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCATCAGTGACCCACATTTGCCG
AGATGTAAATTATGGCTGAATTATCCGATATATACACGCTAACGGAGCTTCTATATTCTTTATATGCCTA
TTCTTACACGTAGGACGAGGCTTATACTATGGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTA
TTTTACTACTCACAGTAATAGCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTG
AGGAGCAACAGTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATTCTACCATTTA

After:
>OK087550.1 Procyon lotor isolate Plo1 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087549.1 Procyon lotor isolate Plo17 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTCGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA
GTAATTACAAATCTCCTGTCAGCTATCCCCTATATCGGATCTAACCTTGTAGAATGAATT
TGAGGAGGGTTTTCAGTAGACAAAGCCACCCTAACACGATTCTTCGCATTCCACTTCATT
CTACCATTTATTATCACAACGCTAGCAATAATTCACCTGCTATTCCTACACGAAACAGGA
TCCAATAATCCTGCTGGAATTACATCAGAATCTGACAAAATTCCATTTCACCCATACTAC
ACCATTAAAGACATTCTGGGAATCCTATTCCTTATTTTTATACTTATAGGTTTAGTGCTA
TTTACGCCAGACCTACTAGGAGACCCAGATAACTACACACCCGCCAACCCCTTAAGCACC
CCACCTCACATTAAACCTGAATGATATTTTCTATTCGCCTACGCAATTCTACGTTCCATT
CCCAACAAACTAGGAGGAGTCCTAGCCCTAGTCCTCTCCATCTTAATCCTAATCATCATC
CCACTCCTACATACCTCAAAACAACGAAGCATAATATTTCGGCCACTTAGCCAATGTTCA
TTCTGATTCCTAGTAGCAGACCTCCTCGTCCTAACATGAATTGGAGGTCAACCAGTAGAA
TATCCCTTCATCATCATCGGCCAACTAGCCTCCATCTTCTACTTTATAATCCTCCTGGTC
CT
>OK087548.1 Procyon lotor isolate Plo11 cytochrome b gene, partial cds; mitochondrial;
TCAAACATCTCAGCATGATGAAATTTCGGCTCCCTCCTTGGAATTTGTTTGCTTCTACAG
ATCGCAACAGGTTTATTCTTAGCCATGCACTATACACCAGATACAACCACAGCTTTCTCA
TCAGTGACCCACATTTGCCGAGATGTAAATTATGGCTGAATTATCCGATATATACACGCT
AACGGAGCTTCTATATTCTTTATATGCCTATTCTTACACGTAGGACGAGGCTTATACTAT
GGCTCCTATACATTCTCTGAAACATGAAATATTGGTGTTATTTTACTACTCACAGTAATA
GCAACCGCTTTCATAGGTTACGTCCTCCCATGAGGACAAATATCATTCTGAGGAGCAACA

Clearer than what? Why is using a complex package requiring several lines of code "clearer" than using standard *nix tools and one liners?

@terdon For the same reasons fasta, nexus, &c. are explicit. It is easier to read at a glance.

I find this much harder to read than the standard tools with their simple syntax. I need to grok 6 lines of code here and manually enter the file name twice with no autocompletion. You also don't mention that you need to install biopython, while the other solutions work out of the box on pretty much any *nix system. On the whole, this is longer, more complex, less portable and significantly slower than either the perl or sed one-liners. This took 6.7s compared to 2s for Perl and 1.6s for sed on a ~300MB file. It does work, but let's stick to the technical and remove the opinion :).

@terdon Of course, you are right from this side. But it reminds me the times when I have only 105 bytes of programme memory & try use the address byte of GOTO as hex-command forwarding on it from another GOTO because it was more effective :-)

Heh, fair. But I just honestly don't find this in any way clearer. It certainly took me more seconds to understand this than any of the other answers. You are obviously more used to python than to perl or sed, so of course you will find it clearer. I am familiar with python, perl and sed, and of the three solutions, this was the least clear since it had more text to read. Everyone's experience is different! :)

@terdon As someone who has taught programming to juniors I would bet a significant amount of money that a larger percentage of 'people in IT' will correctly "decode" what this does compared to any regex based solution. The strength of regexes is their flexibility, not their cleanness or readability. Appending a character at the end of a line really isn't playing to the strengths of regexes.

@DavidMulder oh sure, programmers are more likely to be familiar with python, no question about that. Old school bioinformaticians/*nix geeks like myself, less so. In any case, it was the unsubstantiated "clearer" that rubbed me the wrong way, since it is so very subjective and happens to be wrong in my own case.  And it also feels like using a tank to swat a fly: we're loading all sorts of fancy tools that just make this a slower, less efficient solution.

@DavidMulder I'm with @terdon on this one - I don't find this solution clearer at all. It might appear easier to read at a glance, perhaps because record.description += ";" seems rather innocuous, but what happens if we're now presented with a FASTA file without a description field? Well, you need to know that in this case the description attribute is actually assigned the record id/name. So appending a semicolon to the description attribute will actually create a description field in the output file, basically duplicating the record id/name. Probably not what anyone expected or wanted.

Steve's point deserves emphasis. I love Biopython but its implementation of the id versus description attributes on SeqRecord objects never fails to trip me up.

This script loads all sequences into memory, which wouldn't work for huge input.

seqkit replace -p "$" -r ";" my.fasta > my.renamed.fasta

Where -p is the pattern to match the end of the header text $
Where -r is the replacement pattern

https://bioinf.shenwei.me/seqkit/usage/#replace
Seqkit is an extremely powerful toolkit for working with fasta/fastq files. It can a lot with a single cli command.


text:
Active sites of the query protein does not align with the top hit using HHpred
I am having difficulties with an alignment on HHpred and was wondering if anyone had any answer to the following:
I am trying to align my protein (PDB 2VU9, a Botulinum Neurotoxin Type a) with one of the top hits generated from Vast.
I know that both my query protein and the other protein both have Glutamic acid as their active site. But, my protein has its active site on the 234th amino acid, whereas the other protein has the active site on its 230th amino acid.
So my question is:
What would be the reasoning for which the top matches on Vast do not align at the active site with my protein, despite sharing the same amino acid at its active site?
Any response helps!

Abstract
"HHpred is a fast server for remote protein homology detection and structure prediction and is the first to implement pairwise comparison of profile hidden Markov models (HMMs). It allows to search a wide choice of databases, such as the PDB, SCOP, Pfam, SMART, COGs and CDD."
Söding et al. Nucleic Acids Res. 2005 Jul 1; 33(Web Server issue): W244–W248.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160169/

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


text:
How to create drug synergy plot?
I want to replicate the synergy drug response curve with the Python library given in the following paper.
Could you kindly advice on:

Suitable graphics libraries for Python
Code or pseudo code

Apologies because I don't currently have a test code.
I would really appreciate if anyone could guide me in this.
Thank you in advance.
Fig 1a
https://www.nature.com/articles/s41586-022-04437-2/figures/1

The two major graphics packages in Python are matplotlib and seaborn. The very best graphics package (in my opinion - hence its a comment) is R's ggplot2. To do this in Python you usually use either a pandas Dataframe or Numpy. This example looks like curve fitting

for which one of the 2025 combinations ?

https://www.nature.com/articles/s41586-022-04437-2/figures/5 should help ?? To me its just a faith thing

https://pypi.org/project/synergy/  A python package to calculate, analyze, and visualize drug combination synergy and antagonism. Currently supports multiple models of synergy, including MuSyC, Bliss, Loewe, Combination Index, ZIP, Zimmer, BRAID, Schindler, and HSA.

This is a separate question.

@pippo1980 as per understanding between the delta IC50 and Emax of 108,259 combi_id's present in the file shared with figshare

Not sure there is one one way to obtain a Bliss curve (https://www.frontiersin.org/articles/10.3389/fphar.2017.00158/full) wasnt able to find the equation/bliss model they used but I did a superfast speed read of the article

Mmmmh Think I should start from scratch (Ic50 Ec50 single agent ).. then try tonunderstand synergy... this pic through me off https://els-jbs-prod-cdn.jbs.elsevierhealth.com/cms/attachment/43469fda-2391-474f-a97f-4a7307fa6f99/gr1.jpg

@pippo1980 same here :)

here found something in methods:  Conceptually, every point on the Bliss dose response curve is defined as the product between the anchor viability and the corresponding point on the library dose response curve. Shifts in potency (∆IC50) and in efficacy (∆Emax) were calculated as the difference between the observed combination response and expected Bliss (∆IC50 = Bliss IC50 − combination IC50, and ∆Emax = Bliss Emax − combination Emax). ∆IC50 is reported on a log2 scale.

To get you started. Complex editing in matplotlib isn't intuitive.
#!/usr/bin/env python3

from numpy import linspace, sin
# you don't need math its just to get the value of pi for this data
import math
import matplotlib.pyplot as plt

xaxis = linspace(1,1.75*math.pi,10)
drug2 = linspace(1,1.45*math.pi,10)

drug1 = sin(xaxis)
drug2b = sin(drug2)
plt.plot(xaxis, drug1, xaxis, drug2b)
plt.legend(['Drug 1', 'Drug 2'])
plt.show()

The authors look to have curve fitted and for that you would use import scipy.optimize as opt

shouldnt be drug1+drug2 ?

Better written as Expected and Observed but it is entirely hypothetical to demonstrate matplotlib.


text:
How reliable is Alphafold2 for folding random-coil peptide
I have a peptide X of length 50aa. Alphafold2 predicted the structure which takes a helix form and with high average  pLDDT score (>90).
But during actual experiment (e.g. crystallography, spectroscopy), peptide X is determined to have a random coil structure.
My question is how can I resolve the difference between AF2 prediction and crystallography?
To what extend we can use AF2 peptide structure prediction for downstream analysis (e.g. docking, etc)?

Possibility 1: technical error
AlphaFold2 uses MSA to infer residue co-evolution (due to epistasis) and thus proximity. A MSA can become contaminated with incorrect homologues, especially shorter protein. Generally these will have low scores and be filtered out, but if they don't they will come in like a wrecking ball.
Possibility 2: conformational switching
Even though on circular dichroism a random coil and helix could not be more opposites, several proteins adopt both conformations and switch. Amyloids are a classic example of this, switching between a beta sheet and a random coil. Synuclein is another example.
A technically problematic (vide infra) alternative option is that the helical conformation is only found when the protein is bound to another.
Solution: AlphaFold2 testing
Whereas EBI-AF2 provides a single model, the AF2 algorithm can provide multiple possible solutions, which in many cases are actually bona fide alternate conformations. On a protein I tried, the phosphorylations in PhosphoSitePlus (high-throughput) explained wonderfully a second model.
The model with the best score is not the "most normal" conformer: not only such a statement makes no sense, but is the one that best obeys the evolutionary signals.
I would suggest looking at https://github.com/sokrypton/ColabFold
Set the target number of decoys to more than 5 and if all are helices then it's a technical issue. In particular the MSA data would need to be dissected if one really cared. But if nothing comes up... it would be a low reward area of investigation as one would need to discount other less and less likely possibilities.
If the decoys do show alternative folds, i.e. coils, then a few things spring to mind that could be tested in silico (applicable only to some cases):

different pH sidechain protonation
ion binding
phosphorylations
binding to nucleic acids
binding to known protein candidates —not worth exploring if there's no very likely binding partners. If different paralogues bind differ targets it's tricky
binding to membranes

Looking in PyMOL at the location of the high pLDDT scores and after running APBS (Adv. Poisson-Boltzmann Solver) within PyMOL may reveal something in the helical AlphaFold2 model.
Once the conditions which explain the helical model (e.g. bond to DNA as a heterotetramer in the presence of zinc at pH 10 has a very low ∆G), those conditions could be tested in vitro.

thanks for your reply. 1) When you said "A problematic alternative option..." did you mean, it's very unlikely that it happens? 2) How can I tested different pH, ion binding in-silico? 3) For Adv. Poisson-Boltzmann Solver is that an option in AF2?, if not how can I get it?

Technically problematic. Is very common, but protein protein binding requires you to know what the other protein is —string DB etc. will give you loads that don't hold true— and preferably where —XLMS datasets are sparse. ColabFold can do hetero oligomers but paralogues with different targets lose the res' co-evo signal 2. The various tests can be done in Rosetta/PyRosetta but require a lot of coding and know-how. 3. APBS is a plug in PyMOL —I spelt out the acronym just because I'd would have scrambled the letters had I not.


text:
Sequence alignment with MUSCLE
For protein multiple sequence alignment (MSA), I am using MUSCLE (v5.1). The default setting without any specifications gives a 'bad' output even though my proteins are very similar (UniProt sequence and the corresponding PDB SEQRES sequences). 'Bad' is defined as unnecessary gaps are introduced into the sequences.
However when the super5 algorithm is used the MSA output is correct, which is described here https://drive5.com/muscle5/manual/cmd_super5.html
In addition, when I increase the gap penalty this issue can be solved but as far as I know, super4 v5.1 doesn't give you the opportunity to adjust gap penalty.
Documentation describes super5 quite similarly than the default algorithm, it "just" reduces runtime for long sequences.
Can someone please give me a hint what am I missing and/or why super5 outperforms default in terms of accuracy in such "easy" cases?

Hi there, would you mind making your example reproducible? Or at least show us the two alignments? I wonder how it is possible the gap penalty is not adjustable. Is there a reason why you want to use MUSCLE in particular?

In general I would agree with your comment @KamilSJaron . This particular question is an exception  because the underlying algorithm in context is very new and is precisely focused on the observation of the OP. I would not have been aware of the release of Muscle 5 without the OPs question.

Hi!
Thanks for both of you!
Now I get it. To be honest, there's no particular reason for using MUSCLE. I just found it easier to implement it in Python with many cases and the output also seemed easier to parse than lets say BLAST.

HMM Super5 is a hidden Markov model (HMM) where  gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

Explanation To explain that further HMM uses observed outcomes to estimate/optimise the hidden parameters, the gap penalty being one of the two key hidden parameters. Thus the algorithm is exploring a parameter landscape through reiteration generating multi-multiple alignments (if that makes sense) using different gap penalties and considering the error between the alignments (highlighted below). The search strategy is key to the success of the output. Muscle 5 describes 'introducing variations (perturbations) into HMM parameters', I assume this is singularly with reference to super5. The parameter landscape might be quite small because it doesn't seem to be a directional search strategy in any way (there's no memory in it).
The optimisation algorithm, i.e. what is being observed is a metric called the 'pair-wise error (PWE) and pair-wise dispersion (PWD)', so what appears to be happening is the algorithm is attempting to minimise the number of gaps in an alignment and this would be compatible with your results.
Speed Super5 is breaking the alignment up into smaller alignment and solving each 'centroid' before bringing together the 'centroids'. This is an obvious strategy particularly in the era of multiple processors because it overcomes RAM bottlenecks, we are all doing that in any algorithmic context.
Overall The algorithm was not absolutely clear where PWE and the feedback loop is being deployed whether it was in the generation of centroids or bringing together the centroids and its relationship to the guide tree. There was a large reliance on BLOSSUM62 - which didn't quite make sense in the given context because the background was clear about the different substitution matrices and the impact on gap penalties.
Summary Static gap penalty parameterisation is antithesis of what super5 is about, that is why you cannot use them. The speed is clear its using parallel processing and subalignments - there are no surprises there, thats a common strategy and easy to code these days.
Final points

I am clear the ms has not been accepted and is available here https://www.biorxiv.org/content/10.1101/2021.06.20.449169v1.full#ref-10 . I don't want to comment further on the manuscript for several reasons, which I prefer to avoid explaining.
It is better to wait until the paper is published because it will likely be different.
The algorithm itself will unquestionably be subject to simulation studies by other authors.
Simulation studies were absent in the BioRXIV manuscript.
Muscle historically has been at the forefront of methods to increase sensitivity. If you read the Clustal Omega publication they describe the implementation of historical Muscle sensitivity improvements, in other words they must have been approved over the course of time.

Other authors are likely to investigate whether simulation conditions will result in super5 causing a 'false peak' scenario. That is the Achilles heal in any HMM.
There are certainly situations where super5 is useful in automation, that is without question. In some pipelines there is no option of checking the alignments within the pipeline because there are just to many and it stops the pipeline (human eyes the opposite of a data pipeline).
In practice optical annotation is used to overcome limitations post alignment, i.e. to resolve 'errors' in gap penalties and this is desirable because it prevents a phenomena of 'circular mathematics' which is a caution in any alignment method. The author does discuss this phenomena without calling this directly 'circular mathematics'. Any alignment algorithm can erroneously impose itself on the alignment and the metric of 'error' is extremely important in this regard. Optical analysis breaks circularisation. There is a very, very obvious solution to all this without using HMM at one extreme or 'optical analysis' at the other extreme, but one day at a time.
Example of HMM 'Circularisation' is the Achilles heal of my explanations, because I explain stuff, then reexplain it which more likely increases confusion and intersperse everything with waffle. However, I do follow a HMM style algorithm because next time I explain this it will be more honed and punchier. The hidden parameter being 'clarity' and the observed phenomena being the frequency of blank expressions, which is sought to minimise on each iteration. I do hope my search of the parameter landscape is a bit more directional than that used in Super5.

parameters looking at the documentation
-perturb SEED

It is important to do the same alignment using different random number streams. A given seed should result in the same gap pattern. If there are different gap patterns between different 'SEEDS', there's something suboptimal in the HMM and that is not good at all. If you think the Muscle 5 is implicit perfection, you don't need to check it - but its a trust thing. Its good the author made it explicit, i.e. you can set it.
The other parameters are:
-perm PERM
-consiters N
-refineiters N

The second most important parameter in my view is -perm - this is permutations on the UPGMA guide tree.  If different -perm values give different alignments thats not good. It must work reasonably well because the author will have checked it.

Thanks for the detailed explanation!


text:
How to identify [hyper]thermophilic proteins by BLAST search?
I need to find homologs of certain proteins by BLAST search or similar. But I want filter (select) homologs from organisms that are classified as thermophilic or hyperthermophilic. I am wondering if there are suggestions for me how approach to this. Thank you.

It seems like the main issue is finding a list of species to filter your blast results on. I will therefore address this list of species issue.
There are some resources such as Microbe Directory that contain organismal characteristics similar to what you describe. Their main site appears nonfunctional but their github still works and appears to have the full database. This database includes an extreme_environment field, which appears to be what you are going for.
I feel like PATRIC used to have something similar but I can't locate it now.
For more ideas you might look at the Microbe Directory paper to see if they cite other similar databases.

Warning: Hearsay

As the others have said, filtering by a list of taxonomy ids is the
only one. However, there is one extra consideration if this is for wet
lab applications:  expression yields even when codon optimised may be
affected by phylum choice.
Archaeal protein encoded by codon-optimised genes for E. coli generally perform less well than the equivalent bacterial thermophilic protein. I have not read a study on it, but I have heard its tale over and over and had issues myself. GC-content and the amino acid composition do correlate, but I have heard of good results in E. coli with Geobacillus (high-GC thermophilic firmicute) as with Bacillus (high GC mesophilic firmicute) after all. I don't know if its isoelectric points close to pH 7 or what.
But mid GC-content bacterial thermophiles, such as Aquifex or some clades of Betaproteobacteria, may be a better call than Archaea.
So as you need to filter your blast search by a series of phyla, you may want to be selective if for expression in a different host...

You could be onto something, I agree - but I suspect someone would have already noticed. I had a brief look at this problem and the level of convergent evolution between bacterial families is extraordinary looking at the results of metagenomic studies in thermal environments. @MaximilianPress had the easiest way forward. There's a lot more going on, particularly at the protein level.

I don't know, I do honestly think it's one of those many industry anecdotes: things that are of importance to white biotech but not to the academic publication machinery. I'll grant that they mostly involve the world of fermentors —e.g. E. coli clumps with added thiamine which messes up fermentor OD readings, which if as a reviewer I saw as a study I'd end up rejecting for not being scientific enough :D.

Okay, I've probably misread the aims ... as an evolutionary question its hard.

No no. I guessed as the aims: I speculated the OP was interested in searching for high T_M proteins from an engineering perspective as I've done so myself. An enzyme's K_M is close to physiological concentration of the substrate and K_M generally decreases with decreasing temperature, so it's a classic cheat to get enzymes with lower K_Ms.

Got it .. thanks

I recommend use the local BLAST search. You can form the request for desired taxa in Entrez or through NCBI website & download fasta-file onto your disk.

You know any taxa that are classified as hyperthermophilic except the obvious ones? Doing BLAST is trivial. For me the important thing to find the right databases.

I am scarcely sure that such database exists. There was ExtremeDB but it is inaccessible now. You could try to write a letter to the authors, maybe they changed the domain name.


text:
DimPlot - How to highlight cells with identity colors?
I made this wrapper for DimPlot:
highlight_gene_expression <- function(seurat, gene_counts, colors = NULL, cells = NULL){
  Idents(seurat) <- 'cell_type'
  p <- DimPlot(
    seurat,
    raster = FALSE,
    shuffle = TRUE,
    reduction = "umap",
    pt.size = 0.01,
    cells = cells,
    cols = colors,
  )
  p <- p + guides(color = guide_legend(override.aes = list(size = 4), ncol = 1))

  if (!is.null(cells)) { gene_counts <- gene_counts[cells,] }
  expression_level <- {gene_counts |> rowSums()} / {gene_counts |> rowSums() |> max()}
  p[[1]]$layers[[1]]$aes_params$alpha <- expression_level
  p
}

It takes a Seurat object and a dataframe of gene expression levels.  It returns a UMAP with the transparency (alpha) of each point determined by the gene expression level:
highlight_gene_expression(
  seurat, # a seurat object
  trgd_counts, # A dataframe of gene expression levels.  Not important to understand for this question.
  cells = trgd_cells, # A vector of cell barcodes
  colors = colors # c("Erythrocytes"="#556b2f", "HSPCs"="#483d8b", "Proliferating T Cells"="#3cb371", ...)
)

Taking it one step further, I highlight the desired cells instead of just coloring them:
highlight_gene_expression <- function(seurat, gene_counts, cells = rownames(seurat@meta.data), colors = NULL){
  Idents(seurat) <- 'cell_type'
  p <- DimPlot(
    seurat,
    raster = FALSE,
    shuffle = TRUE,
    reduction = "umap",
    cols.highlight = colors,
    pt.size = 0.01,
    sizes.highlight = 0.01,
    cells.highlight = seurat@meta.data[cells,] %>% split(.$cell_type) |> lapply(rownames),
  )
  p <- p + guides(color = guide_legend(override.aes = list(size = 4), ncol = 1))

  expression_level <- {gene_counts |> rowSums()} / {gene_counts |> rowSums() |> max()}
  p[[1]]$layers[[1]]$aes_params$alpha <- expression_level
  p
}

This produces a similar plot, but with unhighlighted points in gray.  I think it looks quite nice:
# Same params as before
highlight_gene_expression(seurat, trgd_counts, cells = trgd_cells, colors = colors)

The problem is that the colors are in the wrong places.  Dimplot accepts my color vector, but it applies it sequentially when used for highlighting.  This means colors are no longer tied to a specific cluster as they were in the previous example.  A cluster will have a different color now depending on which ones are actually highlighted.
This is not a bug.  Rather, it is expected behavior if you look at the documentation for DimPlot.  The cols param corresponds to an identity class.  cols.highlight does not.  But I want it to.  How can I make cols.highlight explicitly define colors for each cluster?


text:
open source tools for CNV analysis from illumina SNP array
What are the best open source tools to analyse copy number variation using Illumina SNP array data...preferably R based tools, but open to other languages as well.
I'm working with leukemia sample and I have Illumina scan data ..in IDAT/GTC file format and i wish to use tools that will utilise the BAF and LRR metrics.
Microarray platform used- Illumina infinium 850K cytoSNP (beadchip array)

What kind of analysis? What questions are you trying to answer? Should we assume human?

@terdon I am trying to do CNV analysis...to detect copy number gains/loss and copy neutral loss of heterozygosity.... I'm working on acute lymphoblastic leukemia samples

Yes, I mean what kind of analysis? The word "analysis" suggests you aren't just trying to detect CNVs, but need to analyse them. Maybe annotate them? Classify them according to the ACMG criteria? If you just want to detect them, that's a different issue.

Maybe you can try karyoploteR. Great tutorial here: https://bernatgel.github.io/karyoploter_tutorial/


text:
Running into an error when loading Flexaid on PyMol 3.7
I have attached a screenshot from terminal, which is the result of me attempting to run Flexaid on my Mac.
I was wondering if anyone knows why this isn't running?
Thanks

Could you kindly check whether the answer below suffices @microbeprincess?  If so please do mark as "accepted", it helps with site stats

Answer from @maximilian-press, converted from comment. Please improve where possible:
It seems that many of the files are not being found, including both run executables and your pdb files. That's what all the No such file or directory stuff is about.
I would suggest confirming that you have completely installed PyMol and any of the relevant dependencies such as Process_ligand that are not being found. Here is one resource that might be relevant for that tool. Is run_flexaid.sh a script that you have written or a third-party tool? I can't find it by google. Can you show us how you ran it?
I am guessing that sh ./run_flexaid.sh is the terminal command that you are using to run the script, but I still don't know where the script comes from. It's possible that running within the script is not using your environment properly somehow. What happens if you try to run Process_ligand on its own (if that's possible?). Have you confirmed that your input pdb files are in the right directory? Basically, it's hard to know what's going on without seeing the script and knowing where your input files live on your computer.
Note: Comments are temporary "Post-It" notes left on a question or answer. They should be used for clarification or constructive criticism of questions or answers, not for answering (no matter how minor). More details on this StackExchange policy here.


text:
Normalization methods to combine scRNA-seq experiments with different sequencing depths
I am training a classifier to identify a cell type in a particular state of activity using scRNA-seq.  There is a large variation in the sequencing depth (reads average per cell) of the testing data (data I will be classifying not training or validating the classifier with) and my supervisor suggested binarizing the genes.  I was wondering if there are other options I should consider.  I was wondering if there are particular normalization methods for ML classifiers which would be best for this task.

Coincidentally, a recent preprint out of Lior Pacther's lab covers some of these details around scRNA-seq depth normalization.
They find that a "PFlog1pPF" normalization method is able to stabilize variance across genes while still reducing the impact of total cell depth on that variance. If I've read this preprint correctly, "PFlog1pPF" is essentially a transformation of the form $\frac{1}{N_c} \log \left(1 + \frac{k_g}{N_c} \right)$ where $N_c$ is the number of reads for cell $c$ and $k_g$ is the number of reads for feature $g$.
This preprint is mostly dedicated towards differential expression and identifying marker genes in certain cell populations, but that may still help you in your application.
Using a transformation like this (or others that are discussed in this preprint and its references) gives you more granularity and sensitivity for highly expressed genes than what binarizing the data will give you.
And it has additional benefits like stabilized variance, which can cause problems in some downstream applications if you don't keep it under control.
But starting with a binarized input like your supervisor suggested might be a simple first attempt to see if your idea has some potential in the first place.

I don't think you need to complicate the idea of normalisation by introducing machine learning classifiers as a necessary component. Normalisation is common when comparing different datasets for all differential analysis.
Given that you have single cell data, have a look at integration techniques in the Seurat workflows. These include normalisation as a component of integration:
https://satijalab.org/seurat/articles/integration_rpca.html

Its not necessarily about ML, but I am specifically using a classifier, in this case a random forest.  So that may affect one's choice of normalization.  For instance binarization may do just fine for a random forest but you would obviously never use that if you were trying to find DEGs.


text:
Grouping x-axis in scatter plot
I am trying to plot one scatter plot for each group, where each group has continous values in x-axis.
Reference Data
Gives continuous x-axis but all groups are combined.
ggplot(data_for_plotting,
aes(x = log10(Count), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

Gives grouped x-axis but the variation in x-axis is lost.
ggplot(data_for_plotting,
aes(x = factor(Clusters), y = log10(Nonsilent.Mutation.Rate), color=factor(Clusters))) +
geom_point(alpha = 1, size = 1) + theme_bw()

Does achieve the grouping but the gives multiple plots.
ggplot(data_for_plotting,
aes(x = log10(Count),y = log10(Nonsilent.Mutation.Rate),color=factor(Clusters))) +
facet_wrap(~ factor(Clusters), shrink = TRUE) +
geom_point(alpha = 1, size = 1) + theme_bw()

I want x-axis to be grouped based on Cluster where each cluster has continous x-axis values.

You are asking for repeating values in the x-axis. That's the only way you'll maintain a continuous x axis while seeing separate lines for each Cluster. The facet is the right solution - they're not separate plots, they're multiple facets of the same plot.

I'm not sure what you're asking for. The first and the third both seem acceptable from a generic visualization point of view. "Gives multiple plots" seems to be what you want based on the prompt? Ram seems to understand you, maybe it's just me.

@MaximilianPress it's not just you. I'm also having trouble understanding the question.

Possibly draw a cartoon of what you would like to accomplish?

You have two different types of values you want to display on the x-axis: Cluster (a categorical variable) and Count (a continuous variable).
As Ram RS suggested, the easiest way to accomplish this is by taking your first plot (which uses a continuous x-axis scale) and facet along your discrete variable.
This gets you what you want, and the rest of the work is adjusting the theme to make it look how you want.
I'm guessing based on your data that you're trying to create something that looks like this figure from Alexandrov et al., Nature, 2013:

You can accomplish this by customizing the facet and the theme.
Here's an example.
gg <- (
  ggplot(
    data = data_for_plotting,
    mapping = aes(
      x = log10(Count),
      y = log10(Nonsilent.Mutation.Rate),
      color = factor(Clusters)
    )
  )
  + geom_point(alpha = 1, size = 1)
  + facet_wrap(
    ~ Cluster,
    # this ensures that all facets are shown on the same row instead of wrapping
    nrow = 1
  )
  + theme_bw()
  + theme(
    # this removes the spacing between facets
    panel.spacing = unit(0, "lines"),
    # this removes the border around facets
    panel.border = element_blank()
  )
)

There are more parts of the theme you can adjust to your liking, like the facet header (e.g. strip.text).
See ggplot2's reference documentation for details.


text:
Problem of writing scilab code on a project - Morris-Lecar model
This is a computational neuroscience project. We are asked to plot the Morris-Lescar model (question 3 below).
The model should be plotted with scilab or R.
however, I don't know where to start? For the moment, the problem is more in writing the code to represent the equations and systems of equations than in actually solving them.
Any clues as to how to start writing our code would be a huge start!

In concrete terms, we don't know where to start with my working group, even for question 3.a..
We know that nullclines are associated with $dV/dt = 0$ and $dW/dt = 0$, which gives the following entries:
$W = \frac{I-g{Ca}M{\infty}(V)(V-V_{Ca})-g_L(V-V_L)}{g_K(V-VK)}$
$W = W{\infty}(V)$
But here for example, no idea how to write it on scilab, or how to display them graphically by varying I.
I hope the application is now clearer.
Appendix: scilab code.
E_L = -50;
E_Ca = 100;
E_K = -70;
V1 = -1.2;
V2 = 18;
V3 = 2;
V4 = 30;
g_L = 2;
g_Ca = 4.4;
g_K = 8.0;
C = 20;
gamma = 0.04;

minf(v)=.5(1+tanh((v-va)/vb))

winf(v)=.5(1+tanh((v-vc)/vd))

tauw(v)=1/cosh((v-vc)/(2*vd))

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

I tried to make the question more to the point, but I think would be easier to answer if you would be even more specific.


text:
GRanges coverage error with discordant reads
I mapped my fastq files from 4SU-seq experiment to human genome hg19 with bowtie 2 and got this results:
    67716784 reads; of these:
  67716784 (100.00%) were paired; of these:
    29941143 (44.22%) aligned concordantly 0 times
    33244334 (49.09%) aligned concordantly exactly 1 time
    4531307 (6.69%) aligned concordantly >1 times
    ----
    29941143 pairs aligned concordantly 0 times; of these:
      9287621 (31.02%) aligned discordantly 1 time
    ----
    20653522 pairs aligned 0 times concordantly or discordantly; of these:
      41307044 mates make up the pairs; of these:
        32859156 (79.55%) aligned 0 times
        3142531 (7.61%) aligned exactly 1 time
        5305357 (12.84%) aligned >1 times
75.74% overall alignment rate

I also ran a second mapping with --no-discordant and --no-mixed. After that I read the bam files into R like this:
x1 <- BamFile(paste0(bam.path, "14751.OnlyConcordantReads.sorted.bam"), asMates = T)
x2 <- BamFile(paste0(bam.path, "14751.AllReads.sorted.bam"), asMates = T)
x1 <- readGAlignmentPairs(x1)
x2 <- readGAlignmentPairs(x2)
x1_gr <- granges(x1, use.names=T)
x2_gr <- granges(x2,use.names=T, on.discordant.seqnames="drop")

The parameter on.discordant.seqnames is necessary, otherwise an error occurs.
z <- GRanges("chr1",IRanges(236958581,237067281))
xcov <- coverage(x1_gr)
xnum <- as.numeric(xcov$chr1[ranges(z)])
plot(xnum)
ycov <- coverage(x2_gr)
ynum <- as.numeric(ycov$chr1[ranges(z)])
plot(ynum)

I get very unusual coverage values (all summed up to 154 - see image) when I use on.discordant.seqnames="drop". Why is that?

In addition, the coverage does not match my bed files.