text:
Examples of machine learning approaches to validate results where ground truth is lacking?
I am currently looking into some of the published methods for deconvolution of spatial transcriptomics data where each spot does not have single-cell resolution. These methods all rely on cell type annotated scRNA-seq data to build a reference of expected expression for each cell population. So in essence, because ground truth is more or less impossible to obtain for complete cell typing of a spatially intact tissue section, another type of data (scRNA) is used instead. This lets you infer something about the spatial organization of the tissue despite no spatial information being present in the training data, but in the end no ground truth is available for validating the results.
I was trying to think of an analogous machine learning problem, biological or not, where something similar is done to deal with unattainable ground truch. Are there any obvious/good examples of something like this? Maybe some inspiration can be drawn for how such results are validated.

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.

Hi, discussion and opinion-based questions are too generic to usefully fit the Stack Exchange format. Please adjust your question to relate to a specific issue that you're having. If you want to ask general questions for discussion, it would be better to do that on reddit

If the problem concerns spatial clustering of cells then the ML algorithm is without question DBSCAN.
The problem with deep learning and CNN (convolution) not "deconvolution", i.e. the MBE paper cited, is:

The model is an ad hoc construction, which can be also a strength;
Under-training, is resolved by adding more layers;
Overtraining, can only be assessed post-analysis, in ML its algorithm specific;
CNN approaches require large datasets for training of 2 million+.

The issue with ML is its performance tapers at around the 1 or 2 million data set size mark, whereas deep learning can (stress can) continue to use massive data sets to enhance performance (thats its strength).
Whether, DBSCAN has received publications in bioinformatics - I'm not sure. However, what is present in the public domain isn't always the cutting edge work: a lesson I've very recently learnt. DBSCAN is very fashionable for spatial problems at present and in context personally would rank it above CNN.
CNN vs ML The thing is most of us - unless you're working for Google - don't have access to data sets in 2 million+ range. This can be circumvented by transfer learning but it is not sufficiently developed outside Google and a few major players in computer vision. Has Google performed CNN on one million+ cell images - I doubt it. If the problem is ecological species identification via remote cameras - Googles CNN image libraries would be just the ticket (lots of dogs, cats and bear photos) but for microscopy studies I seriously doubt it.
Conclusion
What I am 100% certain on is DBSCAN is VERY fashionable right now and certain ML algorithms repeatedly outperform deep learning in certain situations. What I am trying to say is ML is not necessarily just a low data size version of deep learning: in some instances it surpasses it and I think DBSCAN has that potential. The broader challenges of deep learning hinge around 'transfer learning'. If biology manages to resolve that then deep learning will start surpassing machine learning, but firstly it's about sharing, secondly in e.g. RNN this is not technically trivial.

What you can do is pose a detailed problem in a later question and I could explain how training would be used ... generally in complex analytical problems it is difficult to assess precisely how a solution would be implemented.


text:
Matches, mismatches and indels
I have a bam-file of reads and reference genome. How can I located matches, mismatches and indels in three different files with Phred score in each data of the files?

Try one of these:
GATK:
https://gatk.broadinstitute.org/hc/en-us/articles/360036194592-Getting-started-with-GATK4
freebayes:
https://bioinformaticsworkbook.org/dataAnalysis/VariantCalling/freebayes-dnaseq-workflow.html#gsc.tab=0

source: Stack Exchange

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

sort: votes


text:
Interpretation of differently estimated Polygenic Risk Scores (PRSs) by PLINK2 --score command and population stratification based on these estimates
I would like to point two (out of many) options where PRSs can be estimated using PLINK2 --score command.
Given the Genome Wide Association Studies (GWAS) summary statistics file (generated by PLINK2) is trd.sum and the name given to target binary files in PLINK format (.fam, .bim and .bed) is tmp10k_buu_2 I can estimate PRSs using the following commands.

./plink2 --bfile tmp10k_buu_2 --score trd.sum 3 4 9 header --out trd making use of REF (reference allele)

./plink2 --bfile tmp10k_buu_2 --score trd.sum 3 6 9 header --out trd making use of A1 (alternative allele)

I observe that the PRSs generated by these two options are perfectly negatively correlated, and that is due to the allele used. However, I find a collection of negatives and positives as PRSs of target individuals from (say option 1) and the opposite sign with same magnitude with the other option (say option 2) above.
How to correctly interprete the risk when option 1 and 2 are used separately (say for instance, individual A's PRS estimated from option 1 (PRS_1_A) = -0.552 and individual B's PRS estimated from option 1 (PRS_1_B) = +0.036 and when computed the same using option 2 will give PRS_2_A = +0.552 and PRS_2_B = -0.036. Which individual is at higher risk)? This involves correct population stratification because the target individuals need to be sorted according to the PRS (say low to high risk), then using 2 different options cannot determine the person having lowest risk as per one option gets the highest risk by using the other.
Please clarify correct interpretation of PRSs estimated using the above two options.

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


text:
Input for --score in PLINK 1.9 when generationg PRS using GWAS summary statistics
I'm conducting a simulation and I need to obtain polygenic risk scores (PRS) using genome wide association studies (GWAS) summary statistics.
It is known that the GWAS will give BETA and OR for quantitative and binary phenotypes respectively.
When conducting the simulation for quantitative traits I used the following code to obtain PRS from GWAS performed:
./plink1.9 --bfile binary_file1_name --linear interaction --pheno pheno_file_name --covar cov_file_name --parameters 1 --allow-no-sex
./plink1.9 --bfile binary_file2_name --score plink.assoc.linear 2 4 7 header --out prs_quantitative

This is fine, where plink.assoc.linear file contains BETA.
When conducting simulation for binary traits I used the following code to obtain PRS from GWAS performed:
./plink1.9 --bfile binary_file1_name --logistic interaction --pheno pheno_file_name --covar cov_file_name --parameters 1 --allow-no-sex
./plink1.9 --bfile binary_file2_name --score plink.assoc.logistic 2 4 7 header --out prs_binary

I'm not sure whether we need to transform OR to BETA before inputting the file to --score in performing this for binary traits or, simply use the original plink.assoc.logistic which has OR (instead of BETA) to go ahead.
Simply, do we need to input file containing BETA values "only" for --score or either a file containing OR values or BETA values is fine?

If the GWAS summary statistics file provide SNP effects (BETA), use them directly in PLINK --score to compute PRS. For binary traits, GWAS summary statistics file usually outputs odds ratio (OR), and hence, when constructing PRS using OR, we need to first trasnform OR to BETA (where BETA = log(OR)) and then use the same --score command to generate PRS.
Additional informaion:
Alternatively, this R package (GxEprs) can be used to construct PRSs using GWAS summary statistics for both quantitative and binary traits. If you need to opt out of additional OR to BETA transformation step, this R package is useful, as it implements this transformation within its in-built function PRS_binary() for generating PRSs for binary traits.


text:
RNA folding at specific temperature with ViennaRNA in python
I am trying to get dot-bracket-notations of single stranded RNA via the ViennaRNA python package (https://pypi.org/project/ViennaRNA/) at different temperatures. I have read in the docs (https://www.tbi.univie.ac.at/RNA/ViennaRNA/doc/html/examples_c.html) that the folding prediction at different temperatures is possible with C..but nowhere it is described how to specify the temperature when using python. Does anyone know how to specify the temperature at which the RNA folding is predicted with the ViennaRNA package in python?
Thanks for the help!

I think that the "Global Variables" section of the docs describes this:

For the Python interface(s) SWIG places global variables of the C-library into an additional namespace cvar. For instance, changing the global temperature variable thus becomes

RNA.cvar.temperature = 25

Globals are somewhat surprising to me as a way of handling this. It seems like you should be able to specify this differently. Maybe I'm missing something?

Very cool, thanks.


text:
cromwell files downloaded from a URL have too strict permissions
Cromwell can download files that are URLs (i.e. URLs in the inputs.json).  However, the permissions of these files, once they are downloaded, are very strict (rw------- bmac bmac).  I was wondering if there was a way to control the file permissions of said files; I need them accessible to other users.
Cromwell uses a library called akka.http to download url-files.  As stated in the
cromwell configuration docs:

Cromwell uses akka-http to serve requests. For more advanced
configuration settings, refer to the akka-http documentation.

If I run cromwell with the given files
java -jar ~/cromwell/cromwell-84.jar run fq_count.wdl -i fq_count.url.json
I see the permissions for the url-input file:
ls -l cromwell-executions/fq_count/2a65649a-999e-44ef-bec8-89953e0e3f1f/call-count_reads/inputs/429920399/12673489687364390176
-rw-------  2 bmac  bmac

Question:  Is there a way to control the file permissions, i.e. through akka.http or other?
WDL
version 1.0
workflow fq_count {
    input {
        File fastq_file
    }
    call count_reads {
      input: fastq = fastq_file
    }
}

task count_reads {
  input {
    File fastq
    String memory = "5G"
    String time = "00:10:00"
  }

  command <<<
    set -euo pipefail
    fastq='~{fastq}'
    extension="${fastq##*.}"

    if [ "${extension}" = 'gz' ]; then
  cat_cmd='zcat'
else
  cat_cmd='cat'
fi
# using awk instead of bc as bc isn't in debian slim images
"${cat_cmd}" "~{fastq}" | wc -l | awk '{ print int($1/4) }' > num_reads.txt
  >>>

  output {
    File num_reads = "num_reads.txt"
  }

  runtime {
    cpu: 1
    memory: memory
    time: time
  }
}

Input json
{
  "fq_count.fastq_file": https://portal.nersc.gov/cfs/m342/jaws/test_data/sample.fastq"
}

Is adding a chmod command to the end an option? Ugly, but it might work.


text:
Can you recommend me some programs or R packages for species distribution modelling?
I am working with a whole genus of land snails, trying to predict its distribution in the coming decades and how it might get affected by climate change. I have already worked with MaxENT and the R package Biomod2 (I have tried out all the important algorithms within it as well). I was wondering if any of you can recommend another option thats reliable.

Anything with the negative binomial distribution will always work.

It might be helpful to know what ecological data you are working from.

Not sure if that's what you're asking but I'm working with presence only occurrence data, climate data from worldclim and I also possess phylogenetic data but I have yet to use them in such work. Also @M__ id like it if you were a little more specific and maybe give some reasoning behind your response.

Spatial species distribution always follow a negative binomial distribution at any point of most other factors e.g. climate. When a species enters a human monoculture type environment - then deviation from the negative binomial can occur within that monoculture. The clumping will be approximated via trigonometry: 'cause the clumping can be approximated via a circle, very clumped - small circle, more dispersed bigger circle.
The parameters of the negative binomial are likely to change in relationship to climate and thats what I'd model if I was doing this work.  The 'clumping' (which is what the negative binomial models) is likely to tightened as global warming increases. In practical terms the snails start drying out so they have to clump in the remaining pockets of moisture (this might be more true of slugs than snails ... at a guess). Alot of invertebrates that cannot or don't really fly - are prone to desiccation.
Personally, I go to empirical stats rather than use a package. However, the packages you are using will follow these two principles. So you can do it from first principles. Its fairly simple ecology its not like phylogenetics that have been through decades of development.
If you've DNA from the snails then the phylogenetic Bayesian spatial models around Beast provide a powerful basis for the investigation. That would be very strongly recommended if you've DNA sequence from a good proportion of the snails.


text:
Is there a way I can combine geological data with population genetics data?
So I have ran SAMOVA and made calculations regarding haplotype diversity and nucleotide diversity of several populations. I have also made raster data regarding geographical data like slope, ruggedness etc. Is there a program/package that allows me to combine these to types of data into a map, so I can study if there is a correlation between these indexes and the diversity?

Update I looked into modern use of Mantel's test and it will not accept categorical variables. I now remember the null hypothesis of Mantel's ... the association is random between the matrices, i.e. no genetic structure. Therefore p < 0.05 excluding degrees of freedom (so in practice its smaller).
The options are:

Fst implementing that is easy but interpreting is complicated because it's about what constitutes genetic structure

this can be done in Arlequin https://www.youtube.com/watch?v=U2fj3pfPUlE as a tutorial
Fstat is the more famous package, but notably old

General linear modelling;
Machine learning for example random forest, SVM etc ...

I would use points 1 and point 3, 2 is old school now. I would use point 1 to datamine the ML.
The key to Fst is that 0.2 comprises population structure for complicated reasons. Implementing ML is not trivial but is the trendiest analysis going. It is not simply implementation its checking and verifying the associations and avoiding overtraining, which would be the major criticism.
The ML would permits the genetic structure to be predicted across broad geographic regions and those predictions are tested and used to reformulate the model. However, you'd need significant expertise in ML (and coding).

The technical word for this is allopatry and if that doesn't mean much, please do read on the subject because that is central to the theory behind the calculations you are performing. There is a HUGE bank of theory on the subject and is the central alternative to Darwinian selection, yep it's that big. So this is not a passing calculation but a fundamental calculation at the heart of the model of speciation of the organisms within your study. What you want to show is your species does not demonstrate allopatry but niche selection. Ultimately this may not be precisely decipherable (depends on the geographic distances).
Very briefly the analysis to perform now is the Mantel's test. This should be within the Arlequin package. Mantel's test will assess if there is a 'cline' in your data, i.e. an association with gene frequency and geographic distance - in an allele frequency package (genetic distance of point mutation otherwise).
Geological data versus geographic distance ... Mantel's test will cover both calculations, however I think that a simple numerical code will be needed for categorical data. I'd need to check that carefully, because I get that wrong the calculation with be wrong. Basically you need to tell the matrix correlation test that this categorical data being tested.
The next level of analysis in allele frequency calculations involves plotting the pairwise (not global) Fst against against geographic distance and geological type - if the association in Mantel's test is not linear. If the association is linear the Mantel's will be positive, i.e. null hypothesis will be accept, p > 0.05.
I'll check the location of the calculations in various packages later (not necessarily today). I need to check the categorical variable for Mantel's - I've never done that calculation so am not sure.
What you are looking for is Mantel's to reject geographic distance and accept geological data type. That would be a decisive and very publishable result. Fst will be the fall-back to shift between the hypotheses.
If you have loads of SNPs there are phylogenetic tree solutions based on modelling point mutation. The next layer of biophylo-geographic analysis involves packages within the Beast framework. I'm sure that would work here because the number of mutations would be low, also this is really just geographic distance rather than niche. If you've loads of SNPs this would be a good analysis to perform however, but for intraspecifies mtDNA particularly involving cox1 I personally wouldn't think it would be informative. It is perfectly doable though.

A positive Mantel's test is a rare result, but generally allopatry is very common so it is a VERY strong null hypothesis for the data. What I am trying to say is you can't just reject allopatry without testing it. ... well you can but not if you want to publish in e.g. Molecular Ecology.
Good luck with the study and again I'll update the post later. It looks really good and is more than interesting.

I forgot to mention in my original post that I did run a mantel's test, and in fact it's results are what prompted me to check geographical data in the first place. For two of my 4 species there was no correlation between genetic and geographical distance (p value>0.05) , so my thought process is that diversity comes from other geological factors.

I did the test in R using the packages ape vegan and geodist, due to not wanting to complicate things even more with arlequin. I am pretty sure the p value needs to be above 0.005 for it to be significant but I might be wrong.

@Nickmofoe I've answered the question, as far as I practically can.


text:
Do we design primers for single stranded DNA with same criteria as we do for double stranded DNA?
I am trying to design primers for a single stranded DNA virus.
But I am having a dilemma as whether we design the primers same way as double stranded DNA?
do we have to design forward and reverse primers both or just a single primer forward or reverse works?

If you've e.g. a parvovirus you design primers as normal, i.e. must be a pair of primers as for double stranded DNA. Albeit there's only one strand present after the first round of PCR there both strands will be present and without the reverse primer the PCR can't work.
If you want evidence of this, you've had a SARS-CoV-2 diagnosis? Well thats a forward and reverse primer PCR starting with ssRNA converted to cDNA, i.e. single copy and then a pair of primers PCR the resulting locus.

Even long before covid, starting from single-stranded DNA is a standard PCR technique, for example gene cloning from extracted cDNA. The principle of the chain reaction does not change, you always need two primers for amplification. Else, you always amplify from the original template over and over again, but this is not exponential and you'd produce an excess of forward (or reverse) strand product which has no use.

So we will just remove both the denaturation steps from PCR conditions, and the rest of the steps remains same, is it right?

No, everything is the same. Exact same. Get local supervision for lab work, this is a bioinformatics community.


text:
The mathematics of kinship (degrees and distributions of relatedness on individual and group levels)
I have been thinking about to what extent a biological entity’s characteristics can be decomposed into discrete sections of the genome that are parsed and mapped one-to-one to a physical structure. It is fascinating in this regard how indistinguishable biology is from computational linguistics and formal systems theory in general - DNA is a formal language; and we are still trying to unravel the grammar. Right?
I think the following question depends on what type of “grammar” DNA has. I don’t know the term in the field, but I have been using the term “holographic” in my own thinking, when a sequence in a code is interpreted twice - by two different grammars, or two different parts or stages of a grammar, contributing meaning to multiple realizations in the output. (“Holographic”, as in, a 2D surface appears to be 3D so long as the information is encoded on it in such a way to be parsible / reconstructible - the information is not lost, in the mix, but is able to be separated and mathematically resolved (some sort of non-associativity under some function, is what I suspect).
My question is about to what extent “kinship” is indicated by “having similar genes”.
Genes are virtual (I think) - pure information. It seems like the same genes can come into being without them having replicated from the same ancestor (convergent evolution). The presence of the same gene in two biological entities does not strictly imply closer “kinship” when compared to two biological entities without that shared gene. Conversely, we need not expect entities of close kinship to share the same genes (a sibling can inherit one allele of some gene, another sibling the other one, dominant or recessive, etc.)
Perhaps there is an analogy between language and biology, that I haven’t clearly spelled out, yet.
-> Organisms are to DNA what sentences are to “logical form” (in linguistics).
-> DNA is to mechanisms that build DNA (RNA transcription, protein synthesis, etc.) what logical form is to “universal grammar” (or, the actual system cradling the human aptitude for languages, itself).
-> (The odd part): Biological reproduction is when a “sentence” of language directly, physically transforms (or, contributes a part of itself) into an instance of a “grammar” (a valid grammar, as specified by the parameters/constraints/dictates of “universal grammar”). It is intriguing how the analogy gets challenged here, perhaps to be probed further. For now, it could be seen to stretch to how human language can express grammars (like humans can define grammars - constructed languages, as one example - mathematics - or any formal system at all.)
-> Lastly, biological evolution has some analogy to language evolution (perhaps, a fourth “layer” or “grammar” (a specification of means by which something can change - rules, constraints) of the system), in that there is some stochastic generation which gets selected for some (chaotically changing, varying) context (ie, natural selection; as for humans and language, it might be hard to specify the various reasons, including highly social ones, why language features, or cultural features, are either retained or dissipate (Dawkins’s concept of “memes”)).
I think there are two competing dynamics.
One says that similar genes are not “informationally significant” because the “universal grammar” and the transformation rules/rounds/generations on universal grammars themselves (the ways the system can vary from generation to generation, in the space of possible evolutionary trajectories) can easily generate the same structure with multiplicity, independently.
Another competing force says that presence of similar genes indicates a likelihood of kinship because it is improbable that two people who collectively have a high degree of similarity between their total collection of genes, perhaps (this is how ethnicities arise - a group exchanging genes internally, within a single evolutionary path, pass the same genes around amongst each other; the same with human languages).
I wonder if there is any “equilibrium” theory, where the extent to which a system’s propensity to permit random duplicates of a form, is somehow balanced with the (evolutionary) dynamic where memes (genes, language features) cluster with regards to proximity of interaction.
So we at the end of the day have a situation where similar overall gene profiles is actually a probable indicator of “poly”-kinship (a kinship vector where any two entities have some last common ancestor for each of their genes, not just their last common “parent” - and those kinship vectors can have unique distributions, from homogeneous to with wide range - but the same average).
But, there are still outliers where a person is able to have a similar phenotype because of sufficient randomness in the way genes work, without it indicating a high-value (normed) kinship-vector.
(This could apply to linguistics where it would be fallacious to assume Japanese and Turkish must be related to each other just because they share strikingly agglutinative morphology features.)
How is this modeled mathematically, in bioinformatics?
Thank you.

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

The question is too broad and needs to be focussed. Game theory has been applied to population dynamics to explain diverse behaviours of "fitness". Phylogenetics is the tool to answer to your question and has been applied to linguistics.

By game theory I mean understanding equilibrium selection.

Could you please suggest good sub-questions the post could be broken down into? Experienced people in a field often are able to ask better questions because they have more meta-knowledge, but a beginner has many strands of interconnected thought they can’t organize enough yet to separate. I guess my question is about the “formal grammar” of biology (to what extent is it known). It is quite broad. I guess I’ll polish off new sub-questions, one at a time. Yes, that could be a good approach.

Do you have a specific project or analysis goal in mind? Because bioinformatics is such a broad area of study, It's usually better to start off with your actual problem, and find tools or solutions to fix that.


text:
How to find specific types of assemblies for specific species using entrez tools?
How to find specific types of assemblies for specific species using entrez tools?
Task: Trying to specifically find transcriptomes and associated cDNA data for a list of speices.
I can use this oneliner esearch -db assembly -query "{species name }" | esummary to get a list of assemblies and associated data available for any given species.
I can use esearch -db assembly -query "{species name }" | xtract -pattern ENTREZ_DIRECT -element Count to find any specific info nested within the XML tree.
Now I am not sure how to find if the assembly is a transcriptome (which is what I am looking for). I have tried to CRTL +f terms like rna, peptide,amino acids to no avail and I guess there is a chance I am missing something while looking it over manually.
For example here is the output for "zea mays".
Does anyone have any clues? Or should I be using a different db or different tool?

I strongly recommend switching to SRA https://www.ncbi.nlm.nih.gov/sra
Whether the following would work ...
esearch -db sra -query "{species name }" | esummary

I'm not at all sure because SRA has its own toolkit. Maybe NCBI e-utilities e-summary efetch would work on the SRA database. However, the weblink certainly works for transcriptomics.
NCBI has shifted some of its database to be housed off-site onto GCP and AWS and I have a feeling these may not be accessible to e-utilities, but I could be wrong.
Just to note, what I envisage the question meaning is a complete transcriptomic study for a given species (or given sample) under given experimental conditions.

Comments. I better clarify, SRA will be raw Illumina output, it is not assembled into exome cDNA. e-utils is basically for assembled cDNA, it might give you the information - I really don't think it will download it though. So to answer your questions, yes the data will be there on SRA, whether this is the format you want is a different question.

I am actually looking for RNA-seq data for the whole exome. I will try to see if I can find transcriptions and cDNA using SRA. Hopefully off-housing has not affected the API.

@Sudoh Exome, that would be assumed .... there's a format issue. Post-updated


text:
Transcript vs Primary transcript on phytozome
Could someone help me understand what the difference between transcript and primary transcript on phytozme is? For example, this dataset of A.thaliana has "primary transcript CDS" vs CDS.
Off the top of my head primary transcript represents the initial strand of RNA made from DNA and transcripts represent genes post-processing of the initial strand.
The thing that is confusing me is, why is the transcript file larger than the primary transcript file? The transcript file has about 21k more headers than the "primary transcript". The only explanation I can think of is that due to alternate splicing you could end up with a bunch of isoforms that lend more headers to the transcript file. Is that it?

According to a team lead at Phytozome:

Multi-exon genes can produce different transcripts via process of alternative splicing:  inclusion/exclusion of particular exons in the spliced transcript.Designation of one of these transcripts as “primary” will vary based on the genome under consideration:  the primary might be the one with the highest expression under “normal” conditions, it might be the one whose mutations was first studied, or it might simply be the transcript with the longest CDS.

In all cases, Number_of_transcripts >= Number_of_primary_transcripts                    Number_of_primary_transcripts == Number_of_genes


text:
How to change the position of groups colors in DimPlot of Seurat
I have a Seurat object and plotted the Dimplot for UMAP visualization for 2 variables, as shown in the image below.

Now, the problem is that I want the group by variables such as Non-responder and Responder and anti-CLTA4, anti-CLTA4+PD1, anti-PD1 on the top of the UMAP plot and not on the right side. My desired output would look like the following.

The timely help is highly appreciated, many thanks.

If you plan to export the graph anyway, the easiest way could be to edit it with a vector graphics program. The other option is to save the ggplot object returned by Dimplot() and modify it to suit your needs.

Could you please provide an example of Seurat DimPlot() to ggplot ? @Cloudberry

DimPlot returns a ggplot object, so ggplot functions can be applied to it. You can change the legend position like this:
DimPlot(data) + theme(legend.position = "top")
# or
DimPlot(data) + theme(legend.position = c(.1, .9))

Thank you so much, problem solved. I was trying to solve it from DimPlot() method.


text:
Problem with ClinCNV script
I have some problems with ClinCNV:
My input: Rscript /root/ClinCNV/clinCNV.R --normal /root/samples/final.cov --out /root/samples/ -bed /root/samples/fileGC.bed --folderWithScript /root/ClinCNV/
And there is an output: /root/ClinCNV/clinCNV.R: error: Error in getopt(spec = spec, opt = args) : short flag "b" is invalid
I've already ask a question about ClinCNV and the main idea was to run script with absolute pathways.

You are using -bed. Multi-letter flags are usually set using -- and not -. I don't know this particular program but since all other long options use --, and the only option with a b is your -bed, I suspect you just want --bed instead:
Rscript /root/ClinCNV/clinCNV.R --normal /root/samples/final.cov \
   --out /root/samples/ \
   --bed /root/samples/fileGC.bed \
   --folderWithScript /root/ClinCNV/

This is confirmed by the tool's documentation, actually:

--bed - .bed file with coordinates of on-target regions, GC-annotated. Gene anotations in the last column can be used in output.

So yes, you need --bed and not -bed. Because you used -, the tool reads that as a single character flag (-b) with the value ed and doesn't know what to do with it.

On an unrelated note, you seem to be running these things as root which is a very, very bad idea. One bug in the program or one misconfiguration and you can mess up your entire system. I strongly urge you to never run commands as root unless absolutely necessary.

Of course! It's now obvious -b or --bed would work. However -bed will not work @ИльяБетяев. Upvoted.


text:
pyScenic CLI for ctx is giving error: Not a single module loaded
Hi I ran pyScenic's ctx via the command in command prompt:
pyscenic ctx /home/vboxuser/pyScenic/adj.csv \
    /home/vboxuser/pyScenic/mm10__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather \
    --annotations_fname /home/vboxuser/pyScenic/motifs-v9-nr.mgi-m0.001-o0.0.tbl \
    --expression_mtx_fname /home/vboxuser/pyScenic/csk_filtered_scenic.loom \
    --mode dask_multiprocessing \
    --output reg.csv \
    --mask_dropouts \
    --num_workers 20

It gives me a bunch of warnings and then an error, as shown below:
2023-01-17 09:54:52,135 - pyscenic.cli.pyscenic - INFO - Creating modules.

2023-01-17 09:54:52,165 - pyscenic.cli.pyscenic - INFO - Loading expression matrix.

2023-01-17 09:54:55,458 - pyscenic.utils - INFO - Calculating Pearson correlations.

2023-01-17 09:54:55,482 - pyscenic.utils - WARNING - Note on correlation calculation: the default behaviour for calculating the correlations has changed after pySCENIC verion 0.9.16. Previously, the default was to calculate the correlation between a TF and target gene using only cells with non-zero expression values (mask_dropouts=True). The current default is now to use all cells to match the behavior of the R verision of SCENIC. The original settings can be retained by setting 'rho_mask_dropouts=True' in the modules_from_adjacencies function, or '--mask_dropouts' from the CLI.
    Dropout masking is currently set to [True].

2023-01-17 09:54:56,273 - pyscenic.utils - INFO - Creating modules.
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(
/home/vboxuser/mambaforge/lib/python3.10/site-packages/pyscenic/utils.py:244: FutureWarning: Not prepending group keys to the result index of transform-like apply. In the future, the group keys will be included in the index, regardless of whether the applied function returns a like-indexed object.
To preserve the previous behavior, use

    >>> .groupby(..., group_keys=False)

To adopt the future behavior and silence this warning, use

    >>> .groupby(..., group_keys=True)
  df = adjacencies.groupby(by=COLUMN_NAME_TARGET).apply(

2023-01-17 09:54:56,350 - pyscenic.cli.pyscenic - ERROR - Not a single module loaded

I tried uninstalling and reinstalling pyScenic, but it didn't change anything. It shouldn't be pyScenic overall either that is wrong, since I ran pyScenic grnboost command and that worked fine.
It seems like it might be specific to the pyScenic commands. Has anyone run into the same issues?
Thanks!

Background information is here https://pyscenic.readthedocs.io/en/latest/tutorial.html


text:
Is there a tool to visualize RNA secondary structure?
I want to manually fold and pair two sequences (U4 and U6). And I want to make niceish figures.
I haven't found a single program that does that.
thanks

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.

Hi @user3392394 and welcome. The question is difficult to understand. Are you trying to perform two individual RNA structures, or is this two sequences to make a single structure? There are lots of programs for RNA structure, what I suspect is you have a specific project in mind, which is not being clearly communicated. Alternatively are you just looking for a RNA 2D structure viewer?

You could try mxfold2.
https://github.com/mxfold/mxfold2
% mxfold2 predict test.fa
>DS4440
GGAUGGAUGUCUGAGCGGUUGAAAGAGUCGGUCUUGAAAACCGAAGUAUUGAUAGGAAUACCGGGGGUUCGAAUCCCUCUCCAUCCG
(((((((........(((((..((((.....))))...)))))...................(((((.......)))))))))))). (24.8)

The RNA viewer that is particularly good is VA RNA. This stands for "Visualisation Applet for RNA". Its a java app and run it via,
java -jar /Route2Progam/vaRNA/VARNAv3-93-API.jar

Where Route2Progam is the path to where the jar files is stored, all other info in the run code is literal.
Very nice pictures, point and clicky, drag structures into different shapes and easy to colour ... so for a GUI its great.
They recommend PsuedoViewer for psuedoknots, but if that doesn't mean much ... it doesn't matter.


text:
Creating meta cells from scRNA-seq and scATAC-seq for GRN analysis
We have two 10X Genomics scRNA-seq and scATAC-seq datasets generated from similar cells, but in different experiments - separately, not using the multiome platform. Now, we are trying to use Pando (https://github.com/quadbiolab/Pando) to generate a gene regulatory network.
We have successfully converted our CCA-embedded dataset into h5ad format (https://mojaveazure.github.io/seurat-disk/articles/convert-anndata.html), then used bipartite matching (https://github.com/quadbiolab/organoid_regulomes) to receive a .tsv file matching invidiual cells from both datasets - it gives us the barcodes of the cells in both the scRNA and scATAC datasets. Now, we need to actually generate metacells/pseudocells that contain both modalities inside the same cell to be able to use Pando.
According to Fleck et al. 2022 (https://www.nature.com/articles/s41586-022-05279-8), they state

The metacells were integrated using cluster similarity spectrum (CSS)

which would refer to this technique: https://github.com/quadbiolab/simspec/blob/master/vignette/vignette.md
However, that does not make use of the previously generated bipartite-matched cell list.
I'm a little bit lost here - can someone point me in the right direction to how we can generate the metacells containing both modalities, scRNA and scATAC?
Thanks a lot in advance!


text:
Calculate interval of flanking introns
I am looking for a easy way to calculate the intervals of flanking introns based on a transcript annotation file. I have chosen a very complicated method:
bedtools subtract -s -a genes.bed -b onlyExons.bed > introns_hg38.bed
Then I calculated the start and end positions of the transcripts from transcript.gtf annotation file and searched for adjacent intron intervals in R. Unfortunately I don't find all the flanking introns. Perhaps this is due to the 3'UTR or 5'UTR or other sequences that are not introns or exons. Any suggestions?


text:
Which basecaller for nanopore is the best to produce event tables with information about the block size/move table?
I previously used albacore version 2.3.1 to make initial move tables, but then I re-squiggle using Tombo version 1.5.1 to fix the errors.
Example of move table produced by albacore, in which all the moves have the same size:

When it is re-squiggled by Tombo (which produces variable length moves/blocks):

I am starting a new project, and I can repeat the same pipeline, however, my concern is that albacore version 2.3.1 is too old now and I was wondering to ask if there are better pipelines/software for base-calling that produces initial move tables so that I can build on top of.

ONT's newest basecaller, Dorado, can do this for SAM output with the --emit-moves argument. It was added in this commit.
Dorado can be downloaded from here:
https://github.com/nanoporetech/dorado
Here's an example command line:
$ ~/install/dorado/dorado-0.2.1-linux-x64/bin/dorado basecaller --emit-moves \
  ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5 | \
  samtools view -b - > called_DAE_gDNA_reads_all.bam

And a few lines of example output:
@HD     VN:1.6  SO:unknown
@PG     ID:basecaller   PN:dorado       VN:0.2.1+c70423e        CL:dorado basecaller --emit-moves /home/minion/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5
@PG     ID:samtools     PN:samtools     PP:basecaller   VN:1.10 CL:samtools view -b -
@PG     ID:samtools.1   PN:samtools     PP:samtools     VN:1.10 CL:samtools view -h called_DAE_gDNA_reads_all.bam
@RG     ID:460fc479-3b30-4ed5-801f-a82bbe556101_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T10:22:53.256+00:00        PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=460fc479-3b30-4ed5-801f-a82bbe556101 LB:DAE_gDNA     SM:DAE_gDNA
@RG     ID:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T20:28:34.34+00:00 PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=227379b1-c5fc-465d-9e24-eb771133f11e LB:no_sample    SM:no_sample
003d84eb-ce0f-47ea-9c7c-fe939fe8702e    4       *       0       0       *       *       0       110     GTGCCTCCTTGCTTTGGTTCAGTTCTATCCAGCCTTTAGATATCCCAGCCCTCACAGGGTTTGTGTGACCTTGTTGCATTATCATGCGAAACACTTTCGCGTTTTTAGTG  %%%%%(&%%&%%'((''))*.-,-''&&&'&&&(%$$$%&&&'%%%%&(*))'&&&&)-05>6;>>777887++'&&&&&&%%%&*.....0820001:;;;<<87''&&  qs:i:7  du:f:0.526      ns:i:2104       ts:i:10 mx:i:4  ch:i:1037       st:Z:2023-02-09T07:07:00.271+00:00      rn:i:96929      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:98.674     sd:f:22.395     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1
002e48ca-84bf-4419-8370-943d1ca3ac64    4       *       0       0       *       *       0       107     TATGTCTAACTTACTTCGTTCAGTTGTGAAAACTGGGTGTTTAACCTTTGATTCTATCGTGTTTCCCTAGTTTTCGCATTTGACGTGAAACGCTTTCACATTTTGGT     &&'''&%%%%((((*+)***())((')&%%%%$%%%'--,,,,-6996/,+-5ABBBCCFBCE>88:1764424/--,,--&%$$%*08??;;;;<B<:((('()//     qs:i:9  du:f:0.92975    ns:i:3719       ts:i:10 mx:i:4  ch:i:1058       st:Z:2023-02-09T07:07:08.905+00:00      rn:i:62167      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:85.094     sd:f:23.944     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,1,1,0,1,0,1,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,0,1,0,0,1,0,1,0,1,1,0,1,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1,0,1,0,0,1,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,1,0,0
0138d83a-da9f-4468-bafc-b974c4e9eece    4       *       0       0       *       *       0       234     ATGTTTATGTGTTCTTCGTTAGTAAGCCACCGATATTGACTCCAGGGTTTGTGTAACCTTGTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTCGCTTCAAATAGTTTTAATTCAGTTGTCTAAAAATGTTAACATATGGGTATTATTTCTCTAGAACATTTTAGTTGTGTAATTTCTACCATTTAAATGGGGTGAATAAATAACCAAATCATTGAACAACTGGAATG      &&&')'&%$$##$&&''+,-)''&'')**'&%&''))&&&)'((*+8<AEEFGHJIJIJIICA@;843344:254,,,,CEE>====7<<<>>@=?+***'&&''(''&%%%%$%%&)1//..-+,+*'&'(+2+++45<ACCC66677@75556+((((55>@AEFHFHGFEHHJHHFHIKJGFFFKLHZEGEGGIDFDDCDDERKIEGDG>=>5332/'&&&&)**,2242.      qs:i:10 du:f:0.69975    ns:i:2799       ts:i:10 mx:i:4  ch:i:1933       st:Z:2023-02-09T07:07:01.924+00:00      rn:i:77788      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:97.397     sd:f:27.276     sv:Z:quantile   RG:Z:227379b1c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,1,0,1,0,1,0,1,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,1,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,1,0,0,0,1,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,1,0,0,1,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0,0,1,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,1,0,1,1,0,1,0,1,1,0,1,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,1,0,1,1,0,1,0,1,1,1,1,0,1,0,1,1,1,0,0,1,1,0,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,1,1,1,0,0,1,0,1,1,0,1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,1,1,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1

Based on the commit, the move table is attached to the mv tag for each read, indicating an array of 8-bit signed integers (represented as comma-separated values in the SAM output). I'd guess (based on how the calling has worked in the past) that these represent strides comprising a fixed number of consecutive samples, so (as before) it'll need correction to convert the stride width to variable width events.
For the model I'm using here, my guess from the model toml file is that the stride width is 5:
$ grep -B 8 stride ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0/config.toml
[encoder]
type = "serial"
[[encoder.sublayers]]
type = "convolution"
insize = 1
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 1024
bias = true
winlen = 19
stride = 5
```


text:
Coordinate numbering between forward (+) and reverse (-) strand in Tombo
My question is related to the output of Tombo re-squiggling algorithm on nanopore data.
Given the position/coordinate on one strand (+ or -), I want to find the coordinate of its exact location on the opposite strand.
For example, let's say I want to find data/signal on specific coordinates on both strands on Tombo path after re-squiggling. Let's say position 13936 - 14036 on genome on + and - negative strands. I wrote the following code but the sequence obtained from Tombo is not similar at all. I want to know how I can find the coordinate corresponding to a given position on the opposite strand.
# one read from the + strand
mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> start
13936
>>> end
13942
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_pos.loc[rd_id, :]
>>> mapped_start
11205
>>> mapped_end
20084
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> forward_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2731
>>> win_seq = seq[start_idx: start_idx + 100]
>>> forward_win_seq = forward_seq[start_idx: start_idx + 100]
>>> # one other overlapping read on the other strand (-)
...
>>> mapped_start, mapped_end, read_id, fast5_path, _, _ = reads_per_win_neg.loc[rd_id, :]
>>> mapped_start
10979
>>> mapped_end
20091
>>> f = h5py.File(fast5_path, 'r')
>>> tombo_events = list(f['Analyses/RawGenomeCorrected_000/BaseCalled_template/Events'])
>>> reverse_seq = ''.join([e[4].decode() for e in tombo_events])
>>> start_idx = start - mapped_start
>>> start_idx
2957
>>> reverse_win_seq = reverse_seq[start_idx: start_idx + 100]

Now the forward_win_seq and reverse (or even the reverse complement) of it are not similar:
>>> forward_win_seq
'GAGACGGAGCAGACCCATCTGCTACTGCCCTTTCTATAATAACTAAAGTTAGCTGCCCTGGACTATTCACCCCCTAGTCTCAATTTAAAAAGATCCCCAT'
>>> reverse_win_seq
'GGGCCAGGTGCCTGAGATCCATGTTCCATCCTACCTGCCTGACCTGCCCGGCATTGCCAACGACCTCATGTACAGTGCCGACCTGGGCCCCGGCATTGCC'
>>> complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
>>> com_seq = ''.join([complement[base] if base in complement.keys() else base for base in reverse_win_seq])
>>> com_seq[::-1]
'GGCAATGCCGGGGCCCAGGTCGGCACTGTACATGAGGTCGTTGGCAATGCCGGGCAGGTCAGGCAGGTAGGATGGAACATGGATCTCAGGCACCTGGCCC'

When Tombo re-squiggle and aligns a read it adds two attributes to it: mapped_start and mapped_end. Tombo also gives chromosomes and strands.
Imagine you are looking for some value of a signal[s] in position x to y on a chromosome on both strands. In this range, you probably have reads that span the positive strand (forward) and the reads that cover the negative strand (reverse).
Now let a given read that has mapped to a positive strand and covers the entire [x, y], have attributes: pos_mapped_start and pos_mapped_end. Similarly let the reads mapped to negative strands covering [x, y] have attributes: neg_mapped_start and neg_mapped_end.
For reads mapped to + strand:
start_idx = x- pos_mapped_start
end_idx = y - pos_mapped_start + 1
required_signal= read_signal[start_idx: end_idx]

For reads mapped to - strand:
start_idx = neg_mapped_end - y- 1
end_idx = neg_mapped_end - x
required_signal= read_signal[start_idx: end_idx]

Note that x and y are inclusive and the length of the required signal = y - x + 1.


text:
How can I predict the folding of a synthetic protein?
I want to predict the folding of a chimeric antigen receptor sequence. These sequences are completely made out of small parts of different proteins. Do you know if AlphaFold can be used for this?

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.

From @MaximillianPress's response, the answer is no because a chimera has no history the domains evolving together and that is a central part of the algorithms assumptions. Please do remember to mark an answer as "accepted" when the question has been successfully answered.

I'd consult references about AlphaFold. For example, wikipedia has a "Limitations" page on the tool, though the page doesn't mention anything about synthetic proteins.
The hesitation that I would have is that the tool relies to some extent on coevolutionary models of amino acids to infer interactions at both training and prediction steps, so throwing a synthetic protein at it may confuse it. See e.g. slide 10 here.
In an interview, they acknowledge this rather directly:

Now, if you have a protein which is completely new, which does not have an evolutionary history—a synthetic protein, which has just been designed— then our system would not be as accurate as it would be for existing proteins.

I am not totally clear on whether this is due to learning a latent representation of constraints on similar proteins, or whether the tool attempts to directly find homologs. The paper methods suggest to me that it is the former, but I'm not sure. They describe training on only chunks of a big MSA. If that's the case, I'd guess that they still can use information about the natural chunks of your synthetic protein, and aren't expecting to be able to do a full-length protein alignment.
So I think that you can expect the prediction to be relatively low quality. It is still probably much better than any alternative way of predicting structure!
As always, not a bad idea to post an issue on github asking this question.

This is a very good point indeed, it will rely on convolution as part of its algorithm and chimeras will not follow this assumption.


text:
Tree Building Algorithm that treats gaps as deletions
I'm part of a nanopore sequencing experiment that will sequence several generations of viruses. The intent is to perform directed evolution by putting selective pressure on these viruses and tracking the various mutations that occur. See here for the paper this methodology is based on.
I've been working on a software package that can compile consensus sequences from longread nanopore sequencing (see here). However, there's a wrinkle with generating phylogenetic trees from this type of data. It appears that most tree building software ignores indels. This makes sense, as gaps could easily indicate a lack of coverage from illumina sequencing methods or similar issues rather than betray the existence of an insertion/deletion. However, the nature of our sequencing method is such that a gaps really should be evaluated as mutational changes that we want to keep track of in a phylogenetic tree.
We've explored using NextStrain to generate our trees and visuals, but the tree building algorithms they employ too suffer from this gap assumption (see the issue I raised here). These include IQTree, RAxML and PhyML. I have not had much luck finding an algorithm or software that treats gaps as indels when building a phylogenetic tree. Does anyone know of such an algorithm or package I could employ? I figure I could take a stab at writing my own algorithm, but consider that a Plan B resort. No sense re-inventing the wheel and all that.

Various distance measures can do this given an MSA. See for example https://academic.oup.com/bioinformatics/article/28/4/495/212883. As long as you have a transition matrix that deals with gaps it should work fine. You can build trees with e.g. neighbor joining from distance matrices. References on ML methods for dealing with indels can be found in this paper.

What the OP is saying is correct wft mainstream phylogenetics: the reason is they model point mutations - not alignment reliability. There are numerous methods (> 2 contrasting methods) to do what the OP seeks, its quite easy. There's no need to write their own algorithm and would be discourage (account for lots of parameters). Generic NJ is okish (its pairwise), its not quite what the OP wants. The virus cited in the Cell paper is Sindbis virus (alphavirus) ... very interesting. I think I'm probably going to sit out on this question, but I definitely look forward to seeing the paper. Upvoted.

Speaking personally it will take me a good half day to understand the VEGAS system (Cell paper), which - whilst I get the gist: using viruses as replication vectors to accelerate adaptation, in particularly I think it leverages copy number thresholds of capsid proteins; I don't understand the nuance of the wet-lab system.

I need to do a more thorough testing to determine just how effective this answer may be, but I think I've put together a workflow that is close enough to an appropriate answer that I'd post it here.
In my research I came across gappy, (GitHub here), a piece of software that identifies potentially relevant gaps in a phylogenetic analysis and returns a PHYLIP-formatted file that can be run through most tree-building programs. Gappy accounts for indels, but not substitutions, so the output would need to be combined with another analysis to generate the tree I want.
Fortunately, IQTree allows for a partitioned analysis of mixed data (see here), so by feeding the original alignment AND the gappy output through IQTree you should be able to generate a tree that accounts for both substitutions and indels.
Again, more thorough testing is required, but I tested this possibility using an MSA example taken from this website. Saving the FASTA-formatted example as alignment.fasta, and running that file through gappy2, you get an output file called GAPPY2_alignment.fasta_2-inf.phy. Thus if you create a NEXUS-formatted file like the following:
#nexus
begin sets;
    charset part1 = alignment.fasta: *;
    charset part2 = GAPPY2_alignment.fasta_2-inf.phy: *;
end;

and run it through IQTree using the command iqtree -p example.nex, you get the following tree:

+---ENA|CAA23748|CAA23748.1
|
+-------------ENA|CAA24095|CAA24095.1
|
|              +-------------------------------------------ENA|BAA20512|BAA20512.1
+--------------|
               +-----------------------ENA|CAA28435|CAA28435.1

Now, disclaimer: this is the exact same tree you get when you run alignment.fasta through IQTree, hence why more thorough testing is required. The sample size for this example is quite small, and it's therefore not surprising that accounting for indels does not change the outcome in this scenario. However, I think this serves as a solid baseline workflow, where future updates will be more tweaks than complete overhauls, if that makes sense.
Let me know if you see any glaring flaws in my logic, but I think I'm going to close this question. I may provide edits to this answer in the future if significant tweaks/overhauls wind up being done.


text:
ipdSummary command stops working if --identify m5C_TET is included in command, and does not give any m5C info either
I have installed the latest SMRT Tools (ipdSummary version 3.0). My goal is to extract ipdRatio and modified bases. In the ipdSummary --help command :
--identify IDENTIFY   Specific modifications to identify (comma-separated list). Current options are m6A, m4C, m5C_TET. Using --control overrides this option. (default: m6A,m4C)

Accordingly, I am running this command:
ipdSummary sorted_chr21_000.bam --reference $REF_path/$REF_FILE_NAME.fasta --referenceWindow 20:1-200000 --identify m6A,m4C,m5C_TET -j 20 --methylFraction --outfile $save_path/$chrom/$name

My bam file is sorted and relatively small (<200 MB).
But this command never finishes and does not throw any errors. In different stages of working it just stops. I checked different thread numbers and smaller files, but still, the command doesn't work.
However, once I changed 'm5C_TET' to 'm5C' it worked perfectly fine and very fast on the same bam file, but in the output, I could not see any information on 'm5C' modifications. I checked the new command with many other bam files and still the same.
My question is the following:
How can I get at least 'm5C' modifications? Should I change my SMRT Tools version?

I asked the same question on the ipdSummary GitHub page and I got the following response:
"I think that was an experimental option that was never fully working or validated, and it no longer works. We probably need to update the help text. Use the separate tool jasmine if you want to do m5C calling."
jasmine also is a tool in the SMRT toolkit that analyzes the kinetic signatures of cytosine bases in CpG motifs to identify the presence of 5mC.

@Marjan We have Pacbio sequence files and we would like to extract methylation information of these files.I would like to ask if you know any software or are familiar with what steps should be done afterwards.

This would not constitute as an answer, please create a new post with your question. Also, doing some literature reading could give you ideas on how to deal with the data.

You can detect methylation information using ipdSummary with the argument --identify.
The command in my question is detecting methylation too.


text:
IBM SPSS- 2-years survival of cohort
Over the past few months I have been retrospectively analyzing patient files for my master's thesis in medicine. As the median overall survival of my cohort is low, I have decided that I want to analyze the 2-years survival of my population. However, I have not been able to find any literature on how to perform this analysis on SPSS.
Does anyone here have any experience analyzing the median survival over a specific period of time (i.e., 1-year, 2-years or 5-years) using IBM SPSS?
Thank you in advance!
Leon

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.

What you want is Cox regression (proportional hazard) models to perform survival analysis. It is years since I have used SPSS and would encourage consideration for the use of R (Python is a bit more complicated). There is a commercial go-between which is easier to use than R but more computationally efficient than SPSS, but I cannot remember its name at all.
There is a weblink for using Cox regression in SPSS:  https://www.scalestatistics.com/cox-regression.html
I have no idea whether the information is correct, but it looks simple enough to point and click.
Background
A  proportional hazard model is one of numerous survival analysis which model unit risk, e.g. mortality, over time considering the risk e.g. mortality of the baseline population (covariates). Generally a hazard model will comprise a hazard function of the risk of events per unit time against the baseline covariates, such as age, smoking that sort of thing. Thus you wouldn't want to consider the mortality risk of a 25 year old to comparable to that of a 90 year old patient for a given factor, because there's loads of confounders.
At a minimum SPSS could be slow and this is a regression analysis, but you will find out soon enough.

It appears there's a forward prediction being sought here. A simple regression may be able to capture this, but the gist of the question appears to be a precise temporal forward prediction. Precise forward prediction are complicated and whilst there are definitely R packages that do this, at a precision level of forecast e.g. 2 vs 5 year this definitely ain't SPSS, a long way from its capabilities. The OP will need to learn R to do this study, it isn't trivial, or approximate from a hazard style regression. For this style of calculation significant epidemiological and clinical information is required. For example, these are often used in cancer studies, but precisely what the background epi/clinical is will be important as to whether a given R library is relevant.
Summary The answer to the question is SPSS can't do this and learning R is the way forward, because the OP will need to use an established library and this only place I'm aware that such libraries exist. The only alternative is to approximate using classical survival analysis.

Hi there! thank you for your response.
You are absolutely right, it's as easy as that to perform survival analysis. However, this results in a comparison in overall survival, whereas I'm looking specifically for a 2-years survival analysis but I just can not seem to find out how to do this.


text:
How to retrieve Expression Atlas baseline expression results in JSON format?
I'm looking to retrieve Expression Atlas baseline expression results in JSON format, e.g. the results on this page.
Several years ago I did a similar thing but looking at differential expression results, and was able to access the json format for those results here. However, I can't figure out how to access the baseline expression result equivalent.
Does anyone know how to access the baseline expression results in JSON format, or know how I could find out how to do this?

I'd email the support for that group. EBI is pretty good at getting back to you in general: https://www.ebi.ac.uk/support/gxa


text:
different results of RMSD for the same trajectory file
I want to analyze the RMSD of my trajectory files using MdTraj.
I used two different codes for this purpose:
the first code:
rep = 1
traj_name = "step7_%d.dcd" %rep
rmsd_total = []

while (os.path.isfile(traj_name)):
    t = md.load(traj_name, top = 'step5_input.pdb')
    top = t.topology
    protein_sel = top.select('backbone')
    protein_traj = t.atom_slice(protein_sel)
    rmsd = md.rmsd(protein_traj, protein_traj, 0)
    rmsd_without_zero = [i for i in rmsd if i != 0.0]
    rmsd_total.extend(rmsd_without_zero)

    rep+=1
    traj_name = "step7_%d.dcd" %rep

and the second code:
traj_1 = md.load('step7_1.dcd', top='step5_input.pdb')
traj_2 = md.load('step7_2.dcd', top='step5_input.pdb')
traj_3 = md.load('step7_3.dcd', top='step5_input.pdb')
traj_4 = md.load('step7_4.dcd', top='step5_input.pdb')
traj_5 = md.load('step7_5.dcd', top='step5_input.pdb')
traj_6 = md.load('step7_6.dcd', top='step5_input.pdb')
traj_7 = md.load('step7_7.dcd', top='step5_input.pdb')
traj_8 = md.load('step7_8.dcd', top='step5_input.pdb')
traj_9 = md.load('step7_9.dcd', top='step5_input.pdb')
traj_10 = md.load('step7_10.dcd', top='step5_input.pdb')
traj_11 = md.load('step7_11.dcd', top='step5_input.pdb')
traj_12 = md.load('step7_12.dcd', top='step5_input.pdb')
traj_13 = md.load('step7_13.dcd', top='step5_input.pdb')
traj_14 = md.load('step7_14.dcd', top='step5_input.pdb')
traj_15 = md.load('step7_15.dcd', top='step5_input.pdb')
traj_16 = md.load('step7_16.dcd', top='step5_input.pdb')
traj_17 = md.load('step7_17.dcd', top='step5_input.pdb')
traj_18 = md.load('step7_18.dcd', top = 'step5_input.pdb')
traj_19 = md.load('step7_19.dcd', top = 'step5_input.pdb')
traj_20 = md.load('step7_20.dcd', top = 'step5_input.pdb')

combine = traj_1 + traj_2 + traj_3 + traj_4 + traj_5 + traj_6 + traj_7 + traj_8 + traj_9 + traj_10 + traj_11 + traj_12 + traj_13 + traj_14 + traj_15 + traj_16 + traj_17 + traj_18 + traj_19 + traj_20

protein_sel = combine.topology.select('backbone')
protein_traj = combine.atom_slice(protein_sel)
rmsd = md.rmsd(protein_traj, protein_traj, 0)

I believe these two codes should result in the same RMSD. However, I get two different RMSD diagrams out of these two codes?
Can anyone help me what is the problems with my first code? (since I think the second code is more accurate while is more naive)

Hi @amirre this has been authoritatively answered please kindly consider accepting the post and helpful for the other site members.

It's a reference issue. The second term.
Your first block makes a list of numpy arrays which are the RMSD against frame  0 of each of the stacks, while in your second you make a concatenation of frames, so your reference frame 0 is common.
The first code block could be smartened by splitting (abstracting) the file valid name enumeration and the reading/analysis. Namely, generate a list of valid filenames and then analyse. This separation would allow you recycle more easily the code: in the seocnd code block you have the file list problem again, which is instead resolved by a copypaste job —very inflexible. A one liner list comprehension (or, less so, a functools.reduce) is a nice option.


text:
Flux balance Analysis, Result interpretation, solver Solver Heuristics
i am new to modelling of cellular networks, and I have two basic questions, regarding the result interpration.
I am using the COBRAToolbox to do a flux balance analysis of a whole genome scale model of a bacteria. To get things started i replicated the results published in a paper, I have even been in correspondence with the author of the paper which happens to be the person who created the model i am using. So I am 100% sure that all the boundary conditions exchange fluxes and so on are as in the publication. I was able to reproduce this results and they match the experimental data.
The published results specifically mention the flux of a metabolite X of interest and compare the reaction rate of X to literature values and to laboratory experiments.
I did a FVA of X given all the boundaries, the soultion space is between 0 - 15 [mmol/g_dw/h], and the result when doing the FBA is aorund 7 [mmol/g_dw/h]?
(I tried out multiple solvers, showing the same variation)
My questions:

I understood the FVA as a analysis tool which shows me the possible solution space for a given flux , given a system with all limitations and boundaries.

So, how is the solver deciding which of the possible values are displayd, which heuristics are used?

If the scattering for the possible solution space is so huge, what is the quality measure for this results?

Thanks in advance.
Best Regards

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:
Minimap2 detects too many 16S sequences in metagenome
I'd like to extract bacterial sequences from a metagenome where ~99% of contigs are from a host insect. My current protocol minimaps every 16S sequences present in the sample, then blasts all those sequences to create a library of full genomes that I minimap against. I'm using the full metagenome as my reference file. I wasn't sure if that makes a difference, but figured I'd mention it. I ran the following command to retrieve 16S sequences that are present in the sample:
minimap2 All_reads.fastq 16S_reference_library.fasta > 16S_map.paf

Minimap2 retrieves every entry in 16S_reference_library.fasta. It's highly improbable some of those sequences are present, or that 25,000 species are present at all. Are there any parameters I could adjust in Minimap2 that might help with this?
Alternatively, I'm trying to implement Krakenuniq and Centrifuge, but I'm concerned that, because I don't have a reference genome for the host insect, they'll have a very high false positive rate.

Do you mean to be using your reads (fastq file) as the reference rather than the query? Normal syntax is minimap2 reference.fasta query.fastq. What do you mean by "retrieves"? How good are the alignments, are they *? Can you show us a snippet of the PAF?

16S is a very conserved sequence, which is why it's used for targeted phylogenetic analysis; it makes it easy to amplify and analyse. Unfortunately that conservation is an issue with minimap2, which is built around the idea that matching scattered subsequences within a sequence is good enough for identifying matches to the entire sequence. With a highly-conserved sequence, there's a higher chance that subsequences will be identical, meaning that diverged sequences will appear to be the identical from the perspective of the mapper.
Minimap2 can be tweaked to make it less sensitive and more specific, but you'd be better off using a different mapping tool, especially one designed for phylogenetic analysis and conserved sequences. I recommend kraken2.
pre-built 16S indexes for kraken2 can be found here:
https://benlangmead.github.io/aws-indexes/k2
If you really want to map specifically to your own generated metagenome, then you could try LAST (with results filtered through last-map-probs), which I have found to be a little bit better in dealing with specificity issues.
Aside: I notice that you've stated "bacterial sequences from a metagenome". If these are, in fact, metagenomic reads (and not targeted 16S reads), then you would be better off using one of the larger metagenome bacterial databases  (e.g. plusPF or plusPF-8, depending on computer memory) instead of concentrating only on 16S-matching reads. This will provide much better resolution; 16S only really works down to the family level.

You are taking reads as the reference. Given the similarity of 16S, you will find hits for most 16S in other species. You have to filter out poor hits. This is the faster way if you don't have many contigs.
Alternatively, switch the order of reference and reads. But in this case, you need to use a high k-mer occurrence cut off as most of your k-mers will occur >10,000 times. Something like:
minimap2 -c -f30000 -N100 ref-16S.fa reads.fq

This will be fairly slow as there are too many hits in the database but it will be easier to process the output.
You can try both and see what works for you.


text:
Protein-Ligand in PDB
I am a newbie in bioinformatics. In PDB, while searching for proteins, I sometimes get protein with ligand, such as PDB ID 1J1Z. In cases I get pure protein structure, in some cases I get protein-with ligand.
I am having hard time understanding why I can have ligand in PDB structure. I am PDB is supposed to give me pure protein 3D structure. Then how can we have ligand in protein structure??
And sorry if my question is stupid.

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

In some protein crystallography studies, the protein crystals are grown with the substrate/ligand/inhibitor added, thus the resulting structure contains both the ligand and the protein atom coordinates. This enables studying the protein ligand complex to understand binding, specificity, mechanism of catalysis...
On the other hand in some experiments protein crystals are grown without the ligand, thus the resulting structure contains coordinates just for the protein atoms.
On the third hand, the lucky one, your protein, overexpressed in its natural host [sometimes even in different one], co-purify with its own natural ligand, so the resulting structure contains coordinates for both the protein atoms and the sought for ligand.
To understand ligand binding the protein structure without the ligand should be compared to the structure of the protein-ligand complex.
Many published papers report both types of crystals.
For example the crystal structure of pectin methylesterase in complex with hexasaccharide VI: 2NTP is from the publication https://doi.org/10.1038/sj.emboj.7601816 which is the source of six other structures. These structures include the enzyme with different substrates, as well as the active site mutants (D178A) all of which were prepared to understand the molecular basis of substrate specificity. The conclusion after studying all of the mentioned structures is nicely summed in the sentence from the abstract:

The preferential binding of methylated sugar residues upstream of the catalytic site, and demethylated residues downstream, drives the enzyme along the pectin molecule and accounts for the sequential pattern of demethylation produced by both bacterial and plant PMEs.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Is the Bowtie2 MAPQ score dependent on the database index?
I have two very similar genomes in my bowtie2 database (nuc similarity > 90 %). When mapping with bowtie2 against this database, most of the reads that are mapped to one of these genomes have a MAPQ score of 0. Only at the positions that are different between the two genomes I get higher MAPQ scores.
I also created a database with just one of the genomes, inspecting the MAPQ score of a single read. This single read gets a MAPQ score of 24 when the database contains only one genome, but drops to zero when two similar genomes are in the database.
Is this the expected behavior?
My understanding of the MAPQ score computation is, that the alignment score of a read at a specific location is only compared to the score at other positions within the same reference, not across the whole database.
If the alignment scores for a read are then pretty close to each other, it is assumed that the read does not map uniquely, which is reflected in a low MAPQ score.
So my question is, whether the alignment score of a read is compared to other alignments of this read within the same reference OR against alignments with the whole database index?
Changing the database, would then result in different MAPQ scores.

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.

Yes, this is expected.
If a read maps equally well to two different locations in the reference database, then the minimum probability that read maps to another location is 0.5. MAPQ scores are ten times the negative log10 of the probability of a mismatch, so a mismatch probably of 0.5 would lead to a maximum MAPQ score of $-10*log_{10}(0.5) = 3$. If there's any doubt about the mapping (e.g. a mismatch in the sequenced read), or the read could map to more than two locations, then that probability will drop further. I would not be surprised if these MAPQ scores, where the location is impossible to precisely determine, were rounded down to 0.


text:
baseMean threshold
I have an RNA-seq dataset and I am using DEseq2 to find differentially expressed genes between the two groups.  I used pre-filtering to remove any genes that have no counts or only one count across the samples.
furthermore, I want to:

remove genes in low counts by using a base mean threshold.

remove genes that have low counts compared to the rest of the genes.

Question
Is there a common threshold used for the baseMean or a way to work out what this threshold should be?


text:
Pipeline for paired end RNA sequence data to proteins
Forgive the basic question here, but I'm super novice at this ...
I have a series of paired-end RNA fastq files (e.g. SRR10720226_1.fastq.gz and SRR10720226_2.fastq.gz).  I'm trying to figure out what software I need to use to translate those into proteins.
I was trying to use biopython to do this using the translate method.  The problem is biopython doesn't seem to consider both files together, only one fastq at a time.
Which tool should I use to work with the reads together?  Is there anything else I should be doing besides just translating to protein?

This is a much more complicated question than it might seem.
First, you need to understand how RNA-Seq works, and what your data really is. Your "paired-end files" contain reads, which will contain fragments of transcripts.
Since those are only fragments of transcripts, you don't know what frame they are in, so you don't have a single possible translation. And each read would only give you a small part of the protein.
Since the data is for humans, your best course of action would be to first, align the reads to the annotated human genome (or transcriptome). Then, count the number of reads for each gene (or transcript), and decide on a condition to call a gene (or transcript) present in the sample (typically by enrichment between two sets of samples, in your case a threshold might be good enough). Finally, for the genes (or transcripts) present in the sample, find out what protein (or isoform) they correspond to, and if needed retrieve the sequence.
To be clear, you will have to learn a lot more to get that done. For the first step, easier tools might include kallisto or salmon. You would probably want to use biological replicates; SRP237926 seems to have several samples, reading the paper would be necessary to determine whether they are useful.

Thanks for the suggestions.  It's a little unclear to me exactly what kallisto and salmon are doing.  Are they doing alignment like bwa, or something else?  What's the result out of those programs?  It doesn't appear to be a list of proteins.  Sorry again if I'm saying things that don't make sense.  Still trying to get my feet under me here.

Yes, like bwa, they will let you determine what gene a read comes from (after bwa, you would still need to count the reads within genomic loci). But they actually use a totally different approach that gives basically the same result, and are easier to run. I would recommend to read the papers, although this blog post about kallisto gives a nice background.

Any papers you could suggest would be helpful.  Something of a 101 on the subject would be nice.  Seems like there's lots on info out there but it all assumes a certain level of knowledge that I don't have.

For instance, I'm not sure how kallisto fits into the pipeline I'm trying to build.  It seems to output an abundance report and a bam file if using --pseudobam.  But how do I use those and which tools do I use to get to proteins?

If you search "introduction to RNA-Seq", you'll find many articles; I don't have a particular recommendation. You might want to read a few (or watch videos) until you understand how RNA is linked to protein (and DNA), what RNA-Seq is, and how library prep works.

From RNA-Seq, you could in principle, if you had no choice, reconstruct prot sequences. But since RNA-Seq only gives fragments, this is an extremely hard problem. But you don't need to tackle it: since the human genome is annotated, you already know the set of all possible proteins. So all you need to do is find which of the (known) proteins are of interest. Kallisto will allow you to determine the abundance of transcripts in those samples, so you can use this abundance to select the set of transcripts (and proteins) of interest. The actual protein sequences you can download from a database.


text:
How can I get data of single cell RNA sequence with raw count?
I am using dataset GSE85241, but I can't find the read counts of the dataset. It only provides with RPKM values.
How can I find the read counts of a dataset?

Hi @KevisLin I think the answer below is the correct answer, haci is an expert. Please consider accepting the answer

If the authors did not provide the raw counts, and I don't think it is compulsory, unfortunately you will need to download the raw data, align and count.


text:
beginner RNA-seq Replicate papers
have a good R and statistical analysis background (also with machine learning). in addition, I'm a fresh biotechnology grad. I would like to try to replicate some Rna-seq analysis with R papers (with their provided data). Any SHORT (beginner-friendly) papers to recommend?

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

I like Gierlinksi, et al., Bioinformatics, 2015 for this type of problem.
I like this paper because the methods are descriptive and clear, they have lots of figures that you can reproduce along the way, they spend a lot of time discussing "normalization" and quality control, and the data is easy to download.
Even better, the authors spent a lot of time making sure all the samples were split across lanes and flowcells in a consistent way so that both biological and technical differences could be robustly identified.
They also compared a few different methods published at the time for performing differential expression analysis, so you can see how you can process the same data with different statistical models.
Unfortunately, they don't provide the explicit code they run for this paper, so you can't look at their code as a reference to guide you.
But the methods are descriptive enough that I don't think there should be too much difficulty.

I am a beginner too.  I will try the galaxy (training):

https://training.galaxyproject.org/training-material/topics/transcriptomics/tutorials/ref-based/tutorial.html

Here's https://drive.google.com/file/d/1Lwh6ofxi_bb92mB0HoRrMJObE1MmtcMi/view?usp=sharing
It is from the book of Methods of mathematical oncology
https://github.com/ncbi/sra-tools

Suzuki, T., Poignard, C., Chaplain, M., & Quaranta, V. (2021). Methods of mathematical oncology : Fusion of mathematics and biology, osaka, japan, october 26-28, 2020 (1st ed. 2021 ed.). Singapore: Springer Nature Singapore : Imprint: Springer.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
Cromwell - controlling core/mem usage, and CallCaching, while running locally?
I'm sitting on a mountain of microarray (.idat files) data, and have been tasked to analyze it. I chose the MoChA pipeline, developed by the Broad Institute. It runs through Cromwell, of which I know quite little as it is written in Java and I think it's mostly meant to be applied in HPC environments.
Anyway, my research group does not currently have access to any fancy clusters so I'm running it locally. I could get most of the data through after creating a swap file of size 256GB, but for the remaining sets I would need terabytes of RAM... all because the Cromwell engine doesn't limit its resource usage.
Now, there's an option in the Cromwell config file "concurrent-jobs-limit", that would solve this resource problem, but when I attempt to use it, I get errors.

Cromwell version used: 78
Exact command used: https://pastebin.com/70ebMVE5
Entire program log: https://pastebin.com/hQ8U3f86
Configuration file contents (ctrl+F "concurrent-job"): https://pastebin.com/veDikew7

I tried adding the missing key (line 379) but it merely changes the error to https://pastebin.com/dxDHrKdk
I suspect the override of "LocalExample" (line 341) somehow removes invisible default settings... and right about there I'm out of my depth.
I have asked in the Cromwell slack channel already, no responses in two days. I asked in the Broad Institute forums and they said they don't support Cromwell. Sincerely hoping someone here knows enough about Cromwell to help. I'd really really appreciate it. :]
Joel


text:
Understanding protein-peptide docking analysis
I want to explore some potential reactive peptides and seems a common step to conduct some form of docking analysis. I have a question about the theory and a practical question.
Practical question:
An interesting secondary step seems to be to calculate the binding energy with something like PRODIGY, although I'm struggling to understand the .pdb files to run the webserver.
Once I have created a .pdb on either HPEPDOCK or CABS-dock, how do I find the naming convention for the peptides ( Ligand ID ) within the .pdb file? So that I can run in Prodigy-Ligand.
Theory question:
I'm struggling to understand what to take-away from the docking analysis.
I've been using HPEPDOCK ( & CABS-dock, although the former is significantly quicker and seems a good starting place ).
I've been trying to read some papers and see how they use docking analysis. I get the feeling, particularly for peptide-protein ( known receptor ), it is simply a visual inspection to see if it makes sense. Adding the abstract docking score to double check it is low enough ( < -100 ).

I'm struggling to understand further use - most papers seem to simply state they did such analysis and ~ it looked okay, or ~ chose the good ones.
If anyone can recommend some literature that might help me understand, that would be great!

'Reactive peptides' (covalent peptidomimetics) are generally far from a medchemist's first choice.
A ligand is not a peptide, it is a small molecule. A residue in compbiochem is a unit and can be AAs residues, NCAAs residues, ligands or ions. Caffeine has the residue name/id CFF. Therefore, you need to use the first tab: protein-protein.


text:
What is the best approach to do a PCA analysis to distinguish between two groups?
I have two samples A and B and I do have an RNAseq for both of them. Inside each file (for A and B) I do have a set of FPKMs that give me a fold change for each gene (sample vs control). It looks somewhat like that:
gen_name grp01 grp02 grp03 grp04 fold_change
gene1    10    10    5     5     2
gene2    11    11    11    11    1
gene3    5     5     10    10    0.5
gene4    ...

for sample A
and, say
gen_name grp05 grp06 grp07 grp08 fold_change
gene1    100   100   5     5     20
gene2    110   110   110   110   1
gene3    5     5     100   100   0.05
gene4    ...

for sample B
I need to show with PCA that, foe example, gene1 is sample B is much more overexpressed than in sample A, same for gene3, but it looks like gene2 is not very interesting.
What data must be put in PCA analysis?

Using a package is the best approach. Scanpy is popular right now. I do this stuff via first principles, so my advice isn't cool here. Anyway, I'd combine all data into one dataframe, minus log-fold change. I assume log fold change is the opposing sample? The data must have a column to identify A genes from B genes, or denoted within the gene name. PCA does not give gene by gene comparisons, it will be cluster by cluster. Interpreting each clusters is done independently (by mapping the log-fold-change), i.e. what is clustering together and that cluster represents a cohesive expression group.

The alternative is DEG analysis, gene by gene - which you've already done.

Thank you. I want to show that A is different from B based on an expression of gene1 and gene3 (in this test experiment). I want to do this in R but scanpy is a python package. Is there any way to use prcomp() ?


text:
How to extract residues from specific chain of interest from multiple pdb files?
I'm new to biopython & I'm trying to extract residues from specific chain of each file from a list of 316 pdb files. Eg: 1d8t has four chains A, B, C & D of which my interest is to extract just the residues of C chain. As the chain of interest differs for different files, I'm not knowing how to extract specific chain residue from each file and save it separately.
from Bio import *
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Polypeptide import PPBuilder
structure = PDBParser().get_structure('1d8t', '1d8t.pdb')

model = structure[0]
chain = model['C']

for i in chain.get_residues():
    print(i)

With the above code I was able to extract the residues of C chains not sure how to proceed for multiple files in python.

Create a function and run it for each pdb file.
from Bio.PDB import PDBParser

def print_residues(name, pdb_file, chain_key, out_file):
    parser = PDBParser()
    structure = parser.get_structure(name, pdb_file)
    model = structure[0]
    chain = model[chain_key]
    with open(out_file, 'w') as fp:
        for i in chain.get_residues():
            print(i, file=fp)

Then run the function for all pdb files.
print_residues('1d8t', '1d8t.pdb', 'C', '1d8t.txt')
# ...


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

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

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


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

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

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


text:
What is the function of the heteroatoms
I am trying to make a prediction application using the sequence of proteins.
If we think about the 6COU protein, I'm running the code below
pdb = PDBParser().get_structure("6c0u", "6c0u.pdb")

for residue in pdb.get_residues():
    print(residue)

Some of the output looks like below when looking at residues. What is the function of the heteroatoms here, and follow up question;
When reading a paper on Binding Site Detection of Proteins, there is exactly one sentence in the paper.
"We first clean the input structure by removing all heteroatoms and solvent molecules from the protein structure using the Biopython47 library."
What is the meaning of this sentence if you are working on protein sequences, heteroatoms are not already there, what does it mean to delete them?
If you want to check the paper, here

The text you cited mentions "input structure", so they work with structures, not just "protein sequences".
Heteroatoms = atoms marked as HETATM in the PDB and mmCIF formats.
In these formats each atom is listed as either ATOM or HETATM. These are not particularly meaningful categories, but they are there, so they are easy to use. It's a common practice to remove heteroatoms with grep -v ^HETATM.
Note that in the files from the Protein Data Bank not only solvent and ligands are marked as HETATM, but also "non-standard" amino-acids in a protein. An amino acid is considered standard if it can be genetically encoded. So in the files from the PDB selenocysteine is ATOM and selenomethionine is HETATM. Some programs don't follow these rules and output all protein residues as ATOM.
Removing heteroatoms works well for the vast majority of structures that have no modified residues. But in some cases it will remove also some protein residues, which may not intended.
Update: if you are wondering why there are ATOM and HETATM records, I have never came across the rationale. The PDB format
specification has had both types of records since 1976.
Different versions of the format spec use different wording,
but the meaning is always similar:

ATOM records give atomic coordinates for "standard" groups.
HETATM records give atomic coordinates for "non-standard" groups.

(Although standard residues that not in a polymer are also HETATM.)
In the mmCIF format the same record is preserved in _atom_site.group_PDB with the description:

This data item is provided for compatibility with the original Protein Data Bank format, and only for that purpose.

Most of the structures have only standard residues in protein/DNA/RNA polymers, so removing HETATM records removes only ligands and solvent.

What I understood from the paper's goal is that deleting residues where heteroatoms are not empty. The value het appears blank (het=" ") next to 95% residue of protein. When you get rid of these, the protein sequence will shorten. Could they delete the residues here because they sees it as unimportant for predicting pocket (binding) area of protein? I'm also trying to estimate pocket area using protein sequence should I delete it too?

@drorhun In the image in your question the last residue in the protein chain is PHE (resseq=350). The next Residue objects are ligands and waters.


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

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

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

and

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

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