text:
Infer the new (Ban et. al) ribosomal nomenlature (ex. uL53 ) from the ribosomal protein's sequence

I'm trying to programmatically construct a name for each protein subchain in any ribosome from Uniprot in accordance with Ban et. al's 2014 proposal (excerpt given below) using PDB's and Uniprot's APIs.
This is to be a part of a bigger ribosomal analysis suite released next year and we would like to standardize the names as much as possible.
Ban et. al's proposal: [https://bangroup.ethz.ch/research/nomenclature-of-ribosomal-proteins.html]

I'm however not remotely a biologist (I work in compilers) and would really appreciate some clarification on the way nomenclature works across domains of life and how homology plays into it.
This is my naive approach thus far:

I get a ribosome from PDB let's say (3J9M or 5MYJ or any other for that matter) and split it into protein subchains. So on the order of 20-80 proteins.
For each protein subchain, I can get all the information in the world on it from Uniprot: names, sequence, organisms containing it, publications, synonymous names, you name it...
Right now (again, naively) for those proteins that don't already have a new name assigned to them, I take name-synonyms that PDB has plenty of (from various publications, I assume) and scan them with regexes (ex./[LS]\d{1,2}/g) for things like "L8" or "S15".
Thus, for each subchain, I end up with something like "L15" or "S28" in multiple copies (depending on the abundance of synonymous entries in PDB).
The rut i'm in right now is this: having something like "L15" to match it to a new "Ban-nomenclature" name: they provide the nomenclature look-up tables in terms of these three "historical" naming conventions(i.e "human", "yeast", "bacteria") and a taxonomic range while all the I have is the ribosomal structure(ex `` from which I got the subchain. I'm sure there's plenty of information that I can get

Could somebody elucidate to me which assumptions I can and cannot make in this field and how exactly to use the taxonomic range?
Say, (i) can I assume that all archean names are the same as bacterial homologs?
If I draw a eukaryotic cytoplasmic ribosome, do I match its subchains against yeast column or human column (since both are eukaryotic)?

Somebody suggested that I start grouping these subchains based on sequence, not on a domain and I do have access to each protein's sequence, but I am not sure how to go about it. I was thinking to get each of Ban's new names' UniRef90 cluster, let's say, and then start checking which cluster a given subchain belongs to.
But I'm obviously guessing at this point. Any advice or pointers are very appreciated.

Uniprot90 is too limited, you'll likely get better results from  pfam id (domain family annotation) —that is, part of your proteins will be annotated in a group and fusion protein (if any) will have two annotations. Alternatively, (original) COG groupings are really handy for this but aren't available from Uniprot. Then there's always the option of making your own psi-Blast motif on the NCBI website (there's a tactic to it for best results though) and using local blast with it.

@MatteoFerla thanks a ton, this is super informative! A bit too dense, actually.
Rolling with your first suggestion: as i understand it, a "new name" (ex. "uS15") would be represented by a Pfam-group of homologs and then for a given protein with an old-name i would track it to its Pfam group?

Yes. Correct. What you are trying to do is find what genes are in what cluster, which is straightforward. You have given a researched questions as per guidelines, however, there is honestly too much to read for my lazy self (I skimmed and do not care to read the paper, which from the title reminds me of the XKCD about standard), so I am bound to be missing something. But here are my thoughts of 4 pitfalls: (cont'd)

(A) you may have two protein clusters that have a single PFAM, (B) you have gene fusions (one protein, 2 domains of interest), (C) the a names are meant to cover analogues, i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural), (D) you have a very small protein which is poorly conserved, (E) a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).

This is a copy-paste nearly verbatim of comments so that the question has an answer
What you are trying to do is find what genes are in what homology cluster. This is common problem and there are many solution each with some issues.
Uniprot90 is indeed a cluster of homologues but it is too limited. Whereas you require clusters that span all the domains (universal).
COG
The nearly original grouping, is COG (cluster of orthologous genes). It is really handy for this but aren't really available from Uniprot and NCBI (its curator) has been phasing COGs out. In Uniprot records, they are patchily present as its derivative eggNOG (which has varying degrees of "scope", LUCA being the universal one):
<dbReference type="eggNOG" id="COG0451">
<property type="taxonomic scope" value="LUCA"/>
</dbReference>

Pfam
Another option that is well mapped and always present when possible in Uniprot entries is Pfam ids. These are domain family annotations and are also good option. It is a domain fold group. So domains within protein with homologous folds are grouped into one. If a protein has two known domains, it will have two Pfam entries —repeat proteins have loads, but most protein have one domain. In Uniprot XML you will find these are:
<dbReference type="Pfam" id="PF01370">

Similarly there is InterPro, which is more narrow than Pfam, but generally overlaps.
DIY
Then there's always the option of making your own psi-Blast motif on the NCBI website (setting high number of hits, setting the database to refseq or even PDB and doing a few iterations), saving the motif pattern (PSSM) and using local blast with it.
Caveats
In general you may however encounter these pitfalls:

you may have two protein clusters that share a single id
you have gene fusions (one protein, 2 domains of interest). In Pfam the protein will have two ids, with COGs or similar, you'd have only an anamalous length to go by.
you may have wished to find analogues (i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural). No homology clustering scheme can find these.  Say a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).
you have a very small protein which is poorly conserved


text:
Would DESeq2 be appropriate or should I use another tool for differential gene expression?
Hay, this question is sort of an extension of previous post. We have sequenced extra cellular RNA from culture media. Since the number of cells were limited, we did not have enough material for biological or technical replicates, so each sample has only been sequenced once, and there are essentially three groups and each group has 2 subgroups, main groups are control, group A and group B. I think the table below should help -

group
subgroup
Number of unique samples

Control
Day 3
~30

Control
Day 5
~5

Treatment_1
Day 3
~30

Treatment_1
Day 5
~50

Treatment_2
Day 3
~40

Treatment_2
Day 5
~20

We want to compare the gene expression (Or rather abundance of RNA) between samples based on the Day subgroup.
Since the same sample / culture batch has not been sequenced with replicates, and there is variance in number of samples per group, I was wondering if DESeq2 would be a valid option, as the documentation does mention that without replicates, the output from DESeq2 should be considered only for exploratory purpose.
If DESeq2 is not appropriate, which other tool would be better?
EDIT - I came across NOISeq which I will try once I get any output from DESEq2, if anyone has any other tools / papers which can potentially help me out, do let me know!
Thanks!

If I am understanding your question right, you are after Day 3 vs Day 5. In that sense you should be able to test this with a formula like ~ group + number-of-samples + subgroup. This will enable comparing "days" while controlling for the number of samples and treatment. You might also consider adding an interaction term group:subgroup if you think treatment (or lack of it) has differential effects on day 3 and day 5.

We are after same day comparison between groups, sorry if I had not worded it correctly.

When you are lacking replicates, most statistical tools would either complain or even refuse to run the test you would like to perform. Even when a stats tool gives you an output without replicates, you will not be able to generalize your results anyway. In your case, a gene that you get from such an analysis might be the result of biology or just some technical artefact. I don't think your questions can be answered with the data that you have in hand.

I am aware of the caveats you've mentioned, I just want to do the best I can with the data we currently have. I came across NOISeq, but it doesn't seem to be in development any longer, but can work without replicates.

It's very simple. You can fiddle with these unreplicated data and get unreliable results. That takes time and in the end will be unsatisfying. Or you repeat this experiment twice for n=3 and then have a triplicate. That's what I'd do. These rank-based analysis and custom playing with unreplicated data just never yield anything substantial.

Again, I agree with you, but what I'm looking for now is a tool / methodology for rank-based analysis so that we can make the best of the data we already have. I'm looking at NOISeq for now, if you can suggest any tools, please let me know.

Definitely a fan of NOIseq :-) (long story) ... there's criticism it lacks power. What you are looking at is the significance of the observation against Type 1 error. The general criticism of non-parametrics is that they have low power, so few results. Personally, I'd use it and given how widely used DSeq2 is perform the same analysis via Dseq2. I'm not a fan of parametrics in this context, but in this field they dominate and DSeq2 is the best. Thus I'd use a parametric (DSeq2) and non-parametric NOIseq. The justification is you are using one of each. Personally, I do neither, I use machine learning to direct the thresholds of a parametric - but thats outside the public domain.
The time series is more complicated. Its essentially a 2-way problem: within a time point and between time points. Binary log difference is a good analysis between two time points in this situation.
I assume there was treatment between day 3 and day 5 for the treatment samples but not the control.
The issue of replicates. Personally what I would attempt to take a short cut with replicates: if Treatment 1 and 2 have a shared impact on the sample, then treating them as replicates and compare them against the control would be the first stage of analytics. This way if you go to the expense of producing the replicates described here you can be sure you've got a result before conducting further wet-lab work.
Clearly I don't know what Treatment 1 nor Treatment 2 are, but they could be drugs attacking the same metabolic pathway. Basically to do this you need a biological rationale for combining Treatments 1 and 2 to be considered replicates. I've attempted to re-describe this approach to biological analytics below the "-----"
The post-hoc analysis is more complex. The post-hoc analysis is a database operation assessing the relative binary log fold difference between time points for treatment versus control.

@haci likely represent the orthodox approach (e.g. at review level). What I would attempt is to provide biological justification that Treatment 1 and 2 were replicates. That depends on precise justification, i.e. this time series models the impact of generic interruption of this generic class of drug on pathway X, i.e. something that biochemically combines treatment 1 and 2.
The other analysis is a within time point analysis and if you did that between Treatment 1 and Treatment 2 for day 3 and day 5 separately (2 analyses), this would be used as justification if there ain't much different within a time point. In addition if treatment is post day 3, then Treatment 1, 2, and control are replicates, i.e. a triplicate (that's unequivocal). That provides a basis for moving forward.
However, The interaction @haci is describing is likely a formal model because this is a two-way problem. I rarely disagree but you've got to be careful here because its a time-series which aren't  independent events and most "interaction" based models require independence. Thus day 5 is a result of day 3, therefore they are not independent. I've encountered this as an editor and it can go very wrong. Summary analysis between time series - those are independent, analysis between data points within a time-series are not independent and interaction mixes those up.
Note
Time series experiments are always complex and performing wet-lab replicates would keep things simple (but obviously more expensive). So, if you've good analytics support or you're an analyst then you probably do it without further wet-lab, otherwise send it back to the wet-lab and then follow the recipe.

Thanks for your input, I wouldn't mind listening to the NOISeq long story (Once I'm done with my analysis, but you're welcome to PM). The time series does complicate things and inherently Day 3 will be different than Day 5. There are a good number of samples for each group, but lets see what the analysis show.

source: Stack Exchange

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

sort: votes


text:
Plotting PCA from cell counting data, viable way for analysis?
I wanted to help out a colleague and I wanted to know if the things I'm doing are making sense-ish.
The data consists of tissue sections which had been stained for markers. The tissue segments are some small tubules, with around 6 different genes stained. They have counted the number of cells in segment, as well as how many cells are stained with different markers. They have around 30 samples  (Split into three categories, so 10 samples each category) and a ratio of  cells with marker / area of segment for 6 genes. The goal is to check if any marker is prevalent in any sample type.
They are doing their own analysis, and as a fun side project, I thought I will try performing PCA on the generated data. I made PCA plot with the raw values, but couldn't see any clear clustering.
I'm probably missing a normalization step, and at loss on which method would be best. However I fear even if I perform some normalization, would the analysis make sense?
EDIT: The code I used for the analysis is fairly simple,
values # Contains table of ratios
metadata # weather sample treated with chemo or not
# Performing PCA and merging it with metadata for easier plotting
pca_values <- pca(values)
metadata <- merge(metadata, scores(pca_values), by=0)
# plotting with ggplot
ggplot(metadata, aes(PC1, PC2, color=chemo)) + geom_point()

Please show this PCA plot, preferably with cells tagged for the different samples. Showing the script you used to make this PCA plot would also be helpful.

Have you looked at the eigenvalues (variance) for the other PC axes? It could be that the effect you're looking for is small, and only shows up on other axes.
Additionally, you could use stats to see if there is a significant effect. Check out PERMANOVAs, for instance.

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


text:
Plenty of unmapped reads in human RNA-seq data, why does the largest contig from SPADES alignes to human 18sRNA?
Pretty much the title. We have performed RNA-seq using SmartSeq3 and processed it with zUMIs using human genome (GRCh38) as reference and gencode 44 for annotation. We noticed plenty of unmapped reads, and importantly there were plenty of reads in the negative control (RNA-free water). To investigate the matter, I extracted the unmapped reads and used SPADEs (in --rna mode) to assemble them. I got around 3k contigs, and when I blasted the top 3 longest genes (sizing 2.2k, 1.7k and 1.4k in length) in NCBI, the hits for 2.2k and 1.4k sequence was Cutibacterium acnes SZ2 DNA and 1.7k matched with Eukaryotic synthetic construct chromosome 16
Could this indicate? Are the unmapped reads due to bacterial contamination? Technically there should shouldn't be many reads in the water sample, where could they come from?

plenty of unmapped reads => unmapped to what? Human transcriptome? hg19 genomic sequence? To be more certain that your RNA-Seq reads are from other species you may try to use T2T genome assembly and the best annotation you can find to map the original/whole set of your RNA-seq reads.

Unmapped to GRCh38, human transcriptome, with gencode 44. Sorry I missed this detail.

We see in single-cell RNA-seq sometimes up to 50% unmapped reads while data are of pristine quality in all relevant metrics during downstream analysis. Ask yourself whether spending time on investigating this makes any difference. Probably it won't. Maybe it's really some bacterial contamination, or just some random hits. After all, there are billions of bacterial sequences known and with four nucleotides to sample from some random hits are probably expected just by chance. In any case, it won't make an impact if your data that are mapped are fine. I would proceed with downstream analysis and only get back to this if there is any relevant issue downstream.

The technique is supposed to capture only polyadenalated RNA, so curious how the bacterial fragments got in there. I thought the same as you, however in my rush to write the post I forgot to mention that I am seeing this behaviour in negative controls i.e. RNA-free water as well. The water samples have good amount of reads which align to human genome and a decent chunk of unmapped reads. The analysis I performed was from this subset. Should I still be concerned or just carry on with my analysis?

No method is bias-free. Continue and see if mapped data are fine.

Only 25% of the reads are aligning to exon or intronic region, 25% are aligning to intergenic region and 50% are unmapped. Would this number appear concerning to you?

Sounds indeed like a bit much, but again, see whether data are fine downstream and only come back to this if not.


text:
GATK 4.1.4.0 Mutect stats output
I have adapted my pipeline to the new filtering strategies implemented in gatk 4.1.4.0 and while looking at its output I noticed that the stat file generated by mutect2 has a negative number:
callable -1.538687311E9

I did not find the documentation about this file content but in the gatk forum it is mentioned that is should represent the number of bases above a certain coverage (10 by default).
This happened for a WGS sample, while for WES I obtain a more sensible value (well, at least positive and integer, even if maybe a tad high):
callable        8.6884303E7

Has this happened to anyone else?
I compared the resulting calls obtained for the WGS with a previous version of mutect2 and they are comparable-ish (98% of the new "ones" are in common with the "old" ones), but I am still worried that something fishy is happening, expecially considering that the filtering step of the pipeline uses this stats file.

quick followup: seems to be a bug due to an integer overflow (yay Java!), GATK developers have fixed it and I'm trying it out.

This was answered by the OP -- it was a bug that has been fixed in a subsequent release, 4.1.5.0. See the release notes here: https://github.com/broadinstitute/gatk/releases/tag/4.1.5.0

Made callable sites a Long to avoid integer overflow (#6303)

Issue ticket: https://github.com/broadinstitute/gatk/issues/6302

Hey ya, this could be a useful answer if you would provide some details, like the version where the bug is fixed and a link to the changelog. Right now it provides nothing than what OP wrote.

Sure; done. I was just looking to not have it show up as unanswered anymore since the problem was effectively solved. But I get your point.


text:
Is there a aminoacid mutation that forms similar protein?
Is there an amino-acid shift that will 99% of the time end up resulting the same biochemical function/structure and implications? Like Alanine replaced by Leucine will code in 99% of the cases(at any position that first has Alanine) towards proteins that will function the same without big percentual changes in how good it binds or catalyses.
I do have a big amount of mutations. I'd like to filter them out on chance of biological/phenotypical impact, and since some amino acids look like good substitutes by eye (if I look at the molecular structure)

Maybe. For a more detailed answer, please add more details to your question.

What do you mean by "shift", if you mean frameshift the answer is no

A mutation is most often likely to be either neutral (no effect) or destabiling (protein misfolds). Changes in activity require changes at the active site or the entranceway to it —which are often only a handful of residues, so are very uncommon.
So the answer to your question is most mutations keep the activity, although a large fraction break it. While to get a new activity the process is more fluid and is not boolean.
NB. Not all proteins are enzymes, but for simplicity I will concentrate on enzymes.
Neutrality and core vs. surface
The outer residues of a protein, the surface residues, are very tolerant to mutations. Unless they are part of an binding surface (interface) or the change results in a charge change that is. The core residues on the other hand are more sensitive to mutations.
Destabilisation
One can do a scan in silico or in vitro or in vivo to determine the mutational landscape of a protein, which can take various forms.
This being Bioinformatics stack overflow I'll give an in silico example. The following for example is a heatmap I did some time back with Rosetta pmut_scan, which is a classic tool for this and shows in the y axis the residue it was mutated to and on the x axis the position. A larger change in Gibbs free energy of folding) relative to wild type (∆∆G) means it is destabiling, but anything less than 5 kcal/mol is neutral. As you can see many are neutral, but some half are destabiling —with a caveat for proline, which require a backbone change.

In the in vitro setting one speaks of melting temperature (DSC assays) or turnover number over Michaelis constant (enzyme specific assay). In the in vivo setting one speaks of a fitness cost of a mutation, i.e. the protein is less active so the organism grows less.
Drug resistance
Drugs inhibit enzymes by binding the active site. Therefore resistance is often costly. There are some nice papers that look into the fitness cost and mutational landscape of HIV protease or β-lactamase. In this situation a strong epistasis is seen. That is, a first mutation that confers resistance to the drug greatly destabilises the protein, while a second compensates for the first mutation.
Change in function
An enzyme can catalyse a reaction well on a given substrate, because it was evolved this way (technically the encoding gene evolved, but shh). This is the physiological activity. However, enzymes can also catalyse side activities on similar substrates with a lower turnover and higher Michaelis constant. This is especially true for substrates it was not selected against. This is called "enzyme promiscuity". The tolerance of the physiological activity to mutation is called "robustness", while the ability to increase promiscuous activities is "plasticity". Different enzymes have a different balance of the two.
So a mutation may change the strength of the physiological activity, but may result in a new promiscuous activity that can be used as a stepping stone for selection to increase this activity. Therefore, as you can see it is not a binary switch, except for a few very specific case.
Best amino acid for changes?
Proline is a very strange amino acid as it is a secondary amino acid. It lacks a hydrogen on the backbone nitrogen so it cannot form the required hydrogen bond for an α-helix. And can adopt a cis-conformation with minimal penalty (whereas it is about 5 kcal/mol for other chiral amino acids). So changes to this amino acid are often very destabilising.
Protein cores tend to be hydrophobic so changes to the hydrophobic residues are often deleterious, especially when going from a small residue to a larger one.
The remaining residues are frequently found on the surface so may be neutral. However, some are found in the active site and may be able to shift the parameters of the ligand space for that enzyme. A change in charge may result in an increased specificity for a different set of substrates for example. A size change from a larger amino acid to a smaller one may result in a broader range of promiscuous activities —mutations to alanine for example.
So ironically, polar and charged residues as they are frequently found on surfaces are the residues which are most often neutral, but in a small amount of cases they may change activity.


text:
Detect differentially expressing cells
after comments on my original post, I will ask my question again here
I have data (RNA expression values, obtained with multi-channel in situ hybridization) collected from 1mio human cells. For each cell, I have the expression value of a negative control (non-human RNA) and single-cell measures for different RNA species. I know that the expression level of the negative control correlate in a linear fashion with all the other RNA species when there is no real expression (non-expressing cells). For example, when a cell has a high expression value of the negative control, the signal for all other RNA species will also be high. I now want to detect the cells that express an RNA without having a high negative signal, e.g. the "real" expressing cells.
One thing to have in mind is that some RNA species are abundantly expressed in many cells, other are very rare and have a low expression value.
I tried to do this with a linear model (for each of the RNA species, negative control vs. single RNA), computing z-scores, infer p-values and correct for multiple testing. For some RNA species, it works just fine. As soon as I have many cells expressing a certain RNA, the thing gets more tricky because the slope of the linear model doesn't reflect the linear relationship of negative control vs. RNA species in non-expressing cells because it is more shifted towards the "real" expressing cells.
To visualize the problem, I attached two plots. One plot showing high RNA expression many cells (left plot) and one from an RNA species which is expressed on a very low level in few cells (right plot). Blue colored dots are considered as expressing cells according to the linear model. As you can see, for the scenario shown on the left side, many "real" expressing cells (low expression level) are not detected.

Any ideas which statistical model would better fit to detect those cells?
Thank you for your suggestions!

Have you considered plotting the per-cell ratio of RNA species to negative control? What is the downstream goal of all of this?

@DevonRyan I've tried this but the problem is that you still have to manually introduce a signal-to-noise cut-off to define expressing cells, which in my opinion is prone to experimenter bias. The downstream goal is to binary label cells for expression (0 no expression, 1 expression). So rather than doing it manually, I would prefer a model that accounts for the variable expression levels and the variable s/n ratios.

Do you only have a single negative control or do you have more than one? I assume a few of the RNA species (like the second one in your example) are rarely expressed, is that correct?

@DevonRyan I have a single negative control (besides the targets of interest) that is measured in each individual cell. Your assumption is completely correct, some RNA species are rarely expressed and if they are, at a fairly low level (but still detectable and different from the negative control). To note is that both axes are asinh transformed (1=1.17 counts, 3=10 counts).


text:
Construct the Overlap Graph of a Collection of k-mers
I am doing a course on graph algorithms in genome sequencing and the current assignment involves building an overlap graph from a bunch of k-mers. Formally the problem is as follows:
Input: A collection of Patterns of k-mers.
Output: The overlap graph in the form of an adjacency list.
I have a working algorithm that handles most cases except for self-loops. The problem is that I don't know when/why I should be adding self-loops. Any help on this would be greatly appreciated. Below are some examples illustrating my question and my current solution. Thank you!
Example 1:
Input:
ACT
CTT
TTT
Output:
TTT->TTT
ACT->CTT
CTT->TTT
Example 2:
Input:
CCCC
Output:
CCCC->CCCC
Example 3:
Input:
CT
TG
TG
TC
TT
TC
Output:
TC->CT
CT->TC,TG,TT
TT->TC,TG,TT
The goal here is to represent a collection of kmers as an overlap graph where unique kmer will be represented as a node in the graph. Directed edges will connect nodes whose last $k-1$ characters overlap with the other node's first $k-1$ characters. Therefore, repeated kmers that overlap with a other kmer should be reflected as self loops. That's what I understood, however I am looking for a more general rule/solution.
Here is another example to illustrate this.
Input:
CT
TT
TT
TT
TT
TT

Output:
CT->TT
TT->TT

Current code:
/*
 * Creates an overlap graph from a list of k-mers
 *
 * \param patterns      a vector of string k-mers
 *
 * \return a string containing an adjacency list representation of the
 * overlap graph as described in the problem specification
 *
 * Output Format. The overlap graph Overlap(Patterns), in the form of an
 * adjacency list. Specifically, each line of the output represents all
 * edges (u,v) leaving a given node u in the format “u -> v”, where u and v
 * are both k-mers. If a given node u has multiple edges leaving it (e.g. v and w),
 * the destination nodes are comma-separated in any order. For example, if
 * there exist nodes ACG, CGT, and CGG, the resulting line of the adjacency
 * list would be: ACG -> CGT,CGG
 */

#include <algorithm>
#include <iostream>
#include <map>
#include <vector>
#include <set>
#include <string>

std::string overlap_graph(const std::vector<std::string>& patterns) {
    // Overlap Graph of the k-mers collection patterns, in which a vertex exists for
    // each unique k-mer, and an edge (u,v) exists if the last k-1 characters of u are
    // equal to the first k-1 characters of v.

    if (patterns.size() == 1) return patterns[0] + "->" + patterns[0];
    // Build adjancency list O(n^2)
    std::map <std::string, std::vector<std::string>> adj;
    for (const auto& kmer_u : patterns) {
        // Add kmer as key if it doesn't exist in the map and find its neighbours.
        if (adj.find(kmer_u) == adj.end()) {
            adj[kmer_u] = {};
            for (const auto& kmer_v : patterns) {
                if ((kmer_u != kmer_v) &&
                    (kmer_u.substr(1, kmer_u.length() - 1) == kmer_v.substr(0, kmer_v.length() - 1)) &&
                    (std::find(adj[kmer_u].begin(), adj[kmer_u].end(), kmer_v) == adj[kmer_u].end())) {
                    adj[kmer_u].push_back(kmer_v);
                }
            }
        } // Otherwise continue.
    }

    // Convert to string O(n).
    std::string ans;
    for (auto const& key_value : adj) {
        std::string kmer_u{ key_value.first };
        std::vector<std::string> kmers_v{ key_value.second };
        if (!kmers_v.empty()) {
            ans.append(kmer_u + "->" + kmers_v[0]);
            for (std::size_t kmer_v{ 1 }; kmer_v < kmers_v.size(); kmer_v++) {
                ans.append(',' + kmers_v[kmer_v]);
            }
            ans.push_back('\n');
        }
    }
    return ans;
}

void tests() {
    /* Sample Dataset */
    std::vector<std::string> patterns{ {"AAG"}, {"AGA"}, {"ATT"}, {"CTA"}, {"CTC"}, {"GAT"},
                                       {"TAC"}, {"TCT"}, {"TCT"}, {"TTC"} };
    std::string answer{ "TTC->TCT\nCTA->TAC\nAAG->AGA\nCTC->TCT\nTCT->CTA,CTC\nAGA->GAT\nGAT->ATT\nATT->TTC\n" };
    std::cout << answer << std::endl;
    // std::cout << (overlap_graph(patterns) == answer) << std::endl;
    std::cout << overlap_graph(patterns) << std::endl;
}

void solution() {
    std::vector<std::string> patterns;

    std::string pattern;
    while (std::cin >> pattern) {
        patterns.push_back(pattern);
    }
    std::cout << overlap_graph(patterns) << std::endl;
}

int main() {
    //tests();
    solution();
}

I think this is actually a de-bruijn as opposed to an overlap graph. I believe an overlap graph has an edge whenever any positive length suffix of a node overlaps with the prefix of another node.

Also, since you are only concerned with length $k-1$ overlaps, you can find all of the edges in $O(n)$ time, since each node can only have 4 possible outgoing edges. I can extend my answer to include details on this if interested

@Throckmorton Thank you for the feedback - I really appreciate it. However, in the de Bruijn graph, the nodes are the $|k-1|$  unique k-mers whereas in the overlap graph the nodes are the unique $|k|$ k-mers. For instance, the $k=3$ de Bruijn graph for genome AGT will be AG->GT and the overlap graph is AGT->AGT (I am actually having troubles with this part - why did we need to add the self-loop here?). I agree with you that the de Bruijn graph can be solved in $O(n)$  as we don't need to check for duplicates, ie. for TTTT and $k=3$ the graph is TT->TT, TT.

You can define $k$ to be whatever you want. The definition of a de-Bruin graph is a graph just all $k$-mers as nodes and an edge whenever there is a $k-1$ overlap between two nodes. This is exactly what the docstring in your question states.

Also, I don’t think there should be a self edge for AGT, in fact inter example the node “ACT” doesn’t have a self edge either.

And the O(n) time isn’t due to duplicates, it’s because for a single node you can find all out edges in constant time.

You will want to create a self loop for a $k$-mer $s$ whenever $s_1,...,s_{k-1} = s_2,...,s_k$. This will only happen when all of the characters of the $k$-mer are identical, as $s_1=s_2, s_2=s_3,...,s_{k-1}=s_k$.

The usual way to avoid loops in linked lists is to keep a record of the things that have already been used / visited, and stop when an already visited thing has been seen again (or when there is nothing more to visit).


text:
Interpreting MultiQC plot of GATK BaseRecalibration data?
I'm using GATK4 to build a somatic variant calling pipeline. The pipeline uses MultiQC to aggregate quality control data, and one of the QC measures reported is base quality score recalibration from GATK4's BaseRecalibraton. The description in MultiQC says

This plot shows the distribution of base quality scores in each sample before and after base quality score recalibration (BQSR). Applying BQSR should broaden the distribution of base quality scores.

The black line is for my tumor sample and blue line is my normal sample. I don't understand how the before and after distributions are both shown in the plot as I only see 1 distribution per sample. Am I missing something, or should there be a "Post-calibration Count" plot that isn't shown? If nothing is missing, how should I be interpreting the plot?

Did you actually run the ApplyBQSR tool?

I ran both BaseRecalibrator and ApplyBQSR and the BAM headers show ApplyBQSR has been run. Is the post calibration data only stored in the BAM files or is it also in the tables output by BaseRecalibrator? If the post calibration data is only in the BAM files, I'm not sure how MultiQC is supposed to know the data exists...

There may be an option to have GATK output a metrics file.

Turns out two recalibration tables are needed to show pre-recalibration and post-recalibration plots. I previously only had generated the pre-calibration table by running BaseRecalibrator on the BAM that hadn't undergone ApplyBQSR to adjust base scores. Here are the steps to show pre and post recalibration plots in MultiQC:

Run BaseRecalibrator on your initial input BAM. This will create a recalibration table with the original base scores.
Run ApplyBQSR with the recalibration table from step 1 to create a new BAM with the adjusted base quality scores.
Run BaseRecalibrator on the adjusted BAM generated in step 2 to create a new recalibration table with the adjusted base quality scores.
Use both recalibration tables as inputs to MultiQC. This should create both plots.

hey did you method work ? I have done your steps and it does not create both plots. I believe we need to create a csv file report from gatk AnalyzeCovariates. Will try that now

@alex what part does not work? Does MultiQC not create both plots even though you give it both tables as input?

yes correct. the tables were in different directories and didn't work. I then tried to put the tables under the same directory, did not work either. then i noticed both tables had the same 'sample name'. I changed the names of the tables so multiqc would detect 2 'samples' and correctly use the tables as before and after. It also did not work. strangely their old test data for BQSR on github works. things have probaly changed since then and they did not update the code

Not sure what the issue is then. It has been awhile since I used that MultiQC function and eventually dropped those plots from my pipeline. I would ask the creator of MultiQC, Phil Ewels for help on his Gitter. There's instructions for how to contact him Github.


text:
What is a read count?
A very basic question: What exactly is a read count in shotgun sequencing and how is it calculated? I'm particularly interested in understanding the definition of the read count in the context of shotgun metagenomic sequencing though I assume the same concept applies to other types of shotgun sequencing as well. Surprisingly, a Google search did not results in any answers!
In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts". These counts are often normalized to account for genome size or are normalized to convert counts to proportions. Understanding the concept of a read count may help me better understand all these.
Thanks!

To add to the answer from Bioathlete, "reads" in the traditional genomic sense are the string of nucleotide data inferred from the sequencer from a given fragment of genetic material (could be Illumina, Nanopore, PacBio, 10X, etc.).
So read count first applies for all of the "reads" defined above generated by the sequencer. From a NovaSeq, the total number of bases from a single flowcell is up to 3 trillion bases at 150bp each read for the highest throughput flowcell (so a read count of 3,000,000,000,000/150 = 20B total read count).
The read count for a single Nanopore flowcell is roughly (30 billion bases/avg read length), but the read lengths can be SUPER long (100kb compared to 150bp is substantial). Long read technologies are generally used for different analyses than short reads.
So all of the above reads would be defined as "raw reads", meaning they are reads not yet analysed in any way to be used for data analysis. So using illumina as an example, for a single Whole Genome sequencing sample we might have 800,000,000 raw reads. We then use an aligner like bwa mem to map the reads to a single reference genome (or in your case you would use a different aligner to map to multiple genomes). From here we still have the concept of raw reads, but we also now have the total number of mapped reads. A typical good illumina library for a single reference genome will have >90% mapping rate for a human sample.
We can then say we have 800,000,000 raw reads, and 720,000,000 mapped reads with a mapping to a reference genome of 90%. So we can compare "raw read count" to "mapped read count" to get an idea of the quality of our library.
In the case of metagenomics, you would consider your "raw read count" to still be the total number of reads you got off the sequencer, but the "mapped read count" would be the sum of the reads mapped to each individual reference genome to get a total quantity of "mapped reads" which could give you an indication of library quality. I'm not sure what a good metagenomics library "mapped read count" is, but I'm sure you can find that information.
Lastly, we want to get an idea of the depth of coverage for each mapped read to the genome. This will tell us how well our mapped reads actually covered the reference genome or genomes we are interested in.
I hope this was helpful!

Typically read count is the total number of reads going into the analysis. It could be based off single or multiple sequencing libraries.
Also it can be used to describe the number of reads that align to a region of the reference.  Depth or coverage are also terms used in this case.
You will need to be more specific in your use case for an accurate definition.

I modified my question to provide my specific details. I'm not sure if I understand what you mean by "the total number of reads going into the analysis".

"In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts"" .  In this case it is the number of reads assigned to the taxa from the overall input sequencing data.


text:
Making a boxplot or violinplot for several dataframe with different number of samples
I have 3 data frames for treatment of 3 different drugs and in each of them I have the number of mutation types like insertion, deletion, SNP and total of mutations for each patients. In each group I have different number of patients like below
> head(dat1)
                patient  DEL  INS   SNP MNP total
1:    LP6008337-DNA_H06  927  773 40756   0 42456
2:    LP6008334-DNA_D02 1049  799 31009   0 32857
> dim(dat1)
[1] 21  6
>

 head(dat2)
                Patient   DEL  INS   SNP MNP total
1:    LP6008031-DNA_E01 13552 3374 62105   0 79031
2:    LP6005500-DNA_G01   539  500 43451   0 44490
> dim(dat2)
[1] 33  6
>

> head(dat3)
                Patient   DEL   INS    SNP MNP  total
1:    LP6005935-DNA_F03 39168 16739  58095   0 114002
2:    LP6008269-DNA_D08   849   910 103501  11 105271

> dim(dat3)
[1] 106   6

I want to show if the number of each mutational category and total is different for each drug by boxplot, violinplot or any type of visualization
For example this picture is very nice
Can you help me?

To help us help you, please make your data available with dput() next time.
To achieve your goal with ggplot2, you would need all of your data in one data frame and in the "long format".
dat1 <- "
Patient DEL INS SNP MNP total
LP6008337-DNA_H06 927 773 40756 0 42456
LP6008334-DNA_D02 1049 799 31009 0 32857
"

dat2 <- "
Patient DEL INS SNP MNP total
LP6008031-DNA_E01 13552 3374 62105 0 79031
LP6005500-DNA_G01 539 500 43451 0 44490
"

dat3 <- "
Patient DEL INS SNP MNP total
LP6005935-DNA_F03 39168 16739 58095 0 114002
LP6008269-DNA_D08 849 910 103501 11 105271
"

# read data and add a column indicating data frame of origin (drug in this case)
dat1 <- readr::read_delim(dat1, delim = " ")
dat1$drug <- "drug1"
dat2 <- readr::read_delim(dat2, delim = " ")
dat2$drug <- "drug2"
dat3 <- readr::read_delim(dat3, delim = " ")
dat3$drug <- "drug3"

# bind rows
dat <- rbind(dat1, dat2, dat3)

> dat
# A tibble: 6 x 7
  Patient             DEL   INS    SNP   MNP  total drug
  <chr>             <dbl> <dbl>  <dbl> <dbl>  <dbl> <chr>
1 LP6008337-DNA_H06   927   773  40756     0  42456 drug1
2 LP6008334-DNA_D02  1049   799  31009     0  32857 drug1
3 LP6008031-DNA_E01 13552  3374  62105     0  79031 drug2
4 LP6005500-DNA_G01   539   500  43451     0  44490 drug2
5 LP6005935-DNA_F03 39168 16739  58095     0 114002 drug3
6 LP6008269-DNA_D08   849   910 103501    11 105271 drug3

# convert data to long format, ggplot() expects that
dat_long <- tidyr::pivot_longer(dat, cols = 2:6)

> dat_long
# A tibble: 30 x 4
   Patient           drug  name  value
   <chr>             <chr> <chr> <dbl>
 1 LP6008337-DNA_H06 drug1 DEL     927
 2 LP6008337-DNA_H06 drug1 INS     773
 3 LP6008337-DNA_H06 drug1 SNP   40756
 4 LP6008337-DNA_H06 drug1 MNP       0
 5 LP6008337-DNA_H06 drug1 total 42456
 6 LP6008334-DNA_D02 drug1 DEL    1049

ggplot(dat_long,
                aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug)

EDIT: setting scales argument of facet_wrap() to "free" could be helpful
ggplot(dat_long,
       aes(x = name, y = value)) +
  geom_boxplot() +
  facet_wrap(~drug, scales = "free")

Sorry @haci how we show if the groups are significantly different in terms of each feature let's say the number of SNP? I have made one by geom_signif but it compares features internally rather than between groups

@Angel, this should be a separate question but in general people use ANOVA to test the difference of three or more groups. You will have to make sure that your data meet its assumptions though.

Thank you by var.test function I compared them although I am not sure if t test is suitable

If you use a violin plot instead and ggplot2 to make it then:
geom_violin(...not relevant stuff..., scale="count")

will scale each violin plot according to the number of actual entries in it. I often find this nice so small groups look super skinny and people don't focus on random changes in them so much.

Sorry what you mean by small groups look super skinny in my 3 data types?

I didn't say anything about your groups, reread what I wrote.


text:
Any user friendly way to find rare mutations in whole genome raw?
Is there any user friendly way to find rare mutations in the individual human whole genome sequencing raw data? (from Dante, 30x coverage).
To be more specific, I want to find mutations from this paper:
https://docs.google.com/document/d/1EkRMuD6J0-zyMY3MegKhz7Th0hzhzVXg6DWg2Hm6l2Y/edit?usp=sharing
(it's very short, less than one page).
But I'm confused. In their paper they just talk about genes (RAD21, B3GAT2, SMC3, SCN11A , SCN5A, SCN9A, SCN10A, SCN11A, TRPA1), but not mutations? Can we find that diseased-mutations which they talking about in genes, which their list? Or is there not enough data in the paper/study for this?
And if there is not enough data, which data I need to request from authors?
Or those genes have kind a "gold standard sequences" and if that sequence different from standard - there is "diseased" gene?
(I have only very basic genetics and bioinformatics knowledges)

@ATpoint do you have a script set up to systematically find when users cross-post to BioStar or Reddit? (Just wondering how common this trend is)

That isn't a paper, it's just an abstract. You can read the actual paper here: https://doi.org/10.1016/S0016-5085(17)30765-5 but it isn't available for free. Presumably, the paper will give the actual variants found.

You need to run a SNP finding pipeline, filter for SNPs in those genes and then run a SNP annotation pipeline. (You could annotate before you filter but it will take longer).
SNP finding usually has three steps, 1) read alignment to reference, 2) SNP calling and 3) SNP filtering (by quality).
You then bring in gene annotation data to filter SNPs to just those genes.
You then run SNP annotation using tools like VEP, SIFT, PolyPhen, etc.
People use tools like GATK to do this.
Hopefully that's enough information to get you started.


text:
What is the origin of the 'source Gene ID' references given in the 'gene_presence_absence.csv' output of Roary?
I am learning to use Roary for preparing a pan genome for some lactobacillus strains.  In the 'gene_presence_absence.csv' output of Roary (which I view in excel), a 'source Gene ID' is given for each of the strains included in the analysis.  The Roary web page states that this is:

Presence and absence of genes in each sample, with the corresponding source Gene ID.

Is this 'source Gene ID' generated by Roary, or does it relate back directly to the .gff files I used as input?  If I take one of these 'source Gene ID' references and Grep in the directory containing my .gff files I get no matches, so I am assuming it is generated by Roary...?
I had hoped that this reference would link back directly to the gene within the annotation file, but that seems not to be the case?  Should/can I use this reference for anything?

Did you end up getting an answer to this question? I have the same question.

@Sharb - I think i can now answer this question.  I have added an answer.  Let me know if you fee it needs any clarification?

Looking at this again, I have now worked out (and should probably have noticed previously) that the 'source Gene ID' field from the 'gene_presence_absence.csv' file originates from the ID field of the .gff file for the sequence in question.
You can therefore easily tie the results from the 'gene_presence_absence.csv' file back to the annotation from the sequence in question, and therefore to the location of the nucleotide sequence for the gene in question.
So for instance, if you take a 'source Gene ID' from the gene_presence_absence file - say GEICCBEJ_00950 - and use it in a grep command referencing the genome from which it came, grep "GEICCBEJ_00950" my_annotation_file.gff, then you should see the line from the annotation corresponding to this sequence.
The 'source Gene ID' is also referenced in the pan_genome_reference.fa output from ROARY, which gives "single representative nucleotide sequence from each of the clusters in the pan genome (core and accessory)".


text:
How to subset genes and its nested features from a GFF file using a gene list
I would like to subset a GFF file (gene and nested features) from a gene list.
The GFF file looks like this
##gff-version 3
Scaffold_1      JGI     gene    22901   45904   .       +       .       ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24192   24254   .       +       0       ID=PAC4GC:50510902.CDS.3;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510902.CDS.4;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510902.CDS.5;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510902.CDS.6;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510902.CDS.7;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510902.CDS.8;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510902.CDS.9;Parent=PAC4GC:50510902
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510902.CDS.10;Parent=PAC4GC:50510902
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510902.three_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1      JGI     mRNA    22901   45904   .       +       .       ID=PAC4GC:50510903;Name=Genecv11000002m;longest=0;Parent=Genecv11000001m.g
Scaffold_1      JGI     five_prime_UTR  22901   23284   .       +       .       ID=PAC4GC:50510903.five_prime_UTR.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     23285   23423   .       +       0       ID=PAC4GC:50510903.CDS.1;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24031   24062   .       +       2       ID=PAC4GC:50510903.CDS.2;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24198   24254   .       +       0       ID=PAC4GC:50510903.CDS.3;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     24509   24568   .       +       0       ID=PAC4GC:50510903.CDS.4;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37558   37603   .       +       0       ID=PAC4GC:50510903.CDS.5;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37775   37821   .       +       2       ID=PAC4GC:50510903.CDS.6;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     37927   38228   .       +       0       ID=PAC4GC:50510903.CDS.7;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42345   42702   .       +       1       ID=PAC4GC:50510903.CDS.8;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     42798   43343   .       +       0       ID=PAC4GC:50510903.CDS.9;Parent=PAC4GC:50510903
Scaffold_1      JGI     CDS     44798   45079   .       +       0       ID=PAC4GC:50510903.CDS.10;Parent=PAC4GC:50510903
Scaffold_1      JGI     three_prime_UTR 45080   45904   .       +       .       ID=PAC4GC:50510903.three_prime_UTR.1;Parent=PAC4GC:50510903

And a have the target genes in a list such as
Genecv11033552m
Genecv11003131m
Genecv11036683m
Genecv11012576m
Genecv11003654m
Genecv11012587m

I know that is possible to subset the gff using grep -f gene_list.txt <gff_file>. However this extracts only the gene and mRNA features, missing the CDS and UTR entries, while I would like to subset the gene together with all its children features (mRNA, five_prime_UTR, CDS, three_prime_UTR).
This happens because the ID in CDS and UTR features are the same of the mRNA and not as in the gene feature.
Any ideas?

I think this simple script that uses pyranges solves it. Half of the below is just setup to make a minimal reproducible example.
# pip install pyranges
# or
# conda install -c bioconda pyranges

import pandas as pd
from io import StringIO
import pyranges as pr

#########
# setup #
#########

gff_name = "gtfo.gtf"
contents = StringIO("""Scaffold_1   JGI gene    22901   45904   .   +   .   ID=Genecv11000001m.g;Name=Genecv11000001m.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:50510902;Name=Genecv11000001m;longest=1;Parent=Genecv11000001m.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:50510902.five_prime_UTR.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:50510902.CDS.1;Parent=PAC4GC:50510902
Scaffold_1  JGI CDS 24031   24062   .   +   2   ID=PAC4GC:50510902.CDS.2;Parent=PAC4GC:50510902
Scaffold_1  JGI gene    22901   45904   .   +   .   ID=Gene2.g;Name=Gene2.g
Scaffold_1  JGI mRNA    22901   45904   .   +   .   ID=PAC4GC:WHATEVZ;Name=Gene2;longest=1;Parent=Gene2.g
Scaffold_1  JGI five_prime_UTR  22901   23284   .   +   .   ID=PAC4GC:WHATEVZ.five_prime_UTR.1;Parent=PAC4GC:WHATEVZ
Scaffold_1  JGI CDS 23285   23423   .   +   0   ID=PAC4GC:WHATEVZ.CDS.1;Parent=PAC4GC:WHATEVZ""")

pd.read_table(contents, sep="\s+", header=None).to_csv(gff_name, sep="\t", index=False, header=False)

genes_to_keep = set(["Genecv11000001m.g"])

############
# solution #
############

df = pr.read_gff3(gff_name, as_df=True)

def all_descendants(df, genes_to_keep):
    old_len, new_len = -1, len(genes_to_keep)
    while old_len != new_len:
        genes_to_keep.update(df[df.Parent.isin(genes_to_keep)].ID.drop_duplicates())
        old_len, new_len = new_len, len(genes_to_keep)
    return df[df.ID.isin(genes_to_keep)]

desired_gff = all_descendants(df, genes_to_keep)
print(desired_gff)
#    Chromosome Source         Feature  Start    End Score Strand Frame                                ID               Name longest             Parent
# 0  Scaffold_1    JGI            gene  22900  45904     .      +     .                 Genecv11000001m.g  Genecv11000001m.g     NaN                NaN
# 1  Scaffold_1    JGI            mRNA  22900  45904     .      +     .                   PAC4GC:50510902    Genecv11000001m       1  Genecv11000001m.g
# 2  Scaffold_1    JGI  five_prime_UTR  22900  23284     .      +     .  PAC4GC:50510902.five_prime_UTR.1                NaN     NaN    PAC4GC:50510902
# 3  Scaffold_1    JGI             CDS  23284  23423     .      +     0             PAC4GC:50510902.CDS.1                NaN     NaN    PAC4GC:50510902
# 4  Scaffold_1    JGI             CDS  24030  24062     .      +     2             PAC4GC:50510902.CDS.2                NaN     NaN    PAC4GC:50510902

This is super nice, and I believe way simpler than the script I had come up with! Many thanks!

My scipt outsources the parsing of the gtf to a library, which is probably why it looks nicer. Much complexity is hidden behind that read_gtf call XD

Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython.
Dealing with nested data in bash or with some fancy awk script is probably going to trickier to write than it's worth.
You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with mRNA in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.
If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.
Repeat until you reach the end of the GFF file or your list of target genes is empty.

I have actually done it in python using the "flow" you suggested. Thank you!

Answer as a python script for this task specifically

That script looks good! My only suggestion is that, depending on the GFF file you're loading, you may get slowed down when trying to load the entire file as a pandas DataFrame. If you are trying to parse a multi-gigabyte GFF file, going line-by-line keeps the memory usage low. Otherwise, what you've done looks good

Precisely, I have already felt that pain :). I am actually working on a better solution for this. Thanks for the help!


text:
Stacked barplot for single cell analysis
I want to create a stacked bar-graph with different cell cycles for each cell type within each condition. I have uploaded the file for it. Don't know how to go about it as now I have another condition (cell cycle). I was able to do the bar graph comparing cell cycle per condition but now I want to compare cell type and cell cycle in each condition
output of dput(head(all.combined@metadata))
structure(list(orig.ident = c("treated", "treated", "treated",
    "treated", "treated", "treated"), nCount_RNA = c(1892, 307, 1348,
    3699, 4205, 4468), nFeature_RNA = c(960L, 243L, 765L, 1612L,
    1341L, 1644L), percent.mt = c(0.211416490486258, 1.62866449511401,
    4.45103857566766, 4.4065963773993, 0.0713436385255648, 3.87197851387645
    ), RNA_snn_res.0.5 = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
    "13", "14", "15", "16", "17", "18", "19"), class = "factor"), seurat_clusters = structure(c(11L, 11L, 5L, 6L, 11L, 13L), .Label = c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"), S.Score = c(0.476893835992198, -0.0200784617568548, -0.0335915198305002, -0.0247184276246385, 0.010785196602457, 0.0190008903712199), G2M.Score = c(0.204441469200986, 0.173804859670862, -0.0313235510969097, -0.0376796363661889, -0.0559526905696905, -0.0122031631356698), Phase = structure(c(3L, 2L, 1L, 1L, 3L, 3L), .Label = c("G1", "G2M", "S"), class = "factor"), old.ident = structure(c(7L,7L, 1L, 4L, 7L, 9L), .Label = c("Fibroblast", "T cell", "Macrophage", "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell", "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("treated_AAACGCTAGCGGGTTA-1", "treated_AAAGGTAAGTACAGAT-1", "treated_AAAGTGAGTTTGATCG-1", "treated_AAATGGACAAAGTGTA-1",
    "treated_AACAAAGGTCGACTTA-1", "treated_AACAGGGTCCTAGCCT-1"), class = "data.frame")

output of dput(tail(all.combined@metadata))
structure(list(orig.ident = c("untreated", "untreated", "untreated",
"untreated", "untreated", "untreated"), nCount_RNA = c(901, 823,
1184, 1835, 1147, 1407), nFeature_RNA = c(482L, 479L, 649L, 1043L,
604L, 709L), percent.mt = c(1.77580466148724, 2.91616038882138,
4.22297297297297, 3.86920980926431, 2.0052310374891, 4.05117270788913
), RNA_snn_res.0.5 = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
"1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12",
"13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    seurat_clusters = structure(c(7L, 7L, 7L, 14L, 7L, 7L), .Label = c("0",
    "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",
    "12", "13", "14", "15", "16", "17", "18", "19"), class = "factor"),
    S.Score = c(-0.0320858200243315, 0.0304725660342869, 0.0215996091745327,
    0.0384166213301423, 0.144956251122548, -0.0242770509986111
    ), G2M.Score = c(0.0904224391544142, 0.050148242050667, -0.0178041670730754,
    -0.0112596867977946, -0.0519554524339088, -0.0136533184257381
    ), Phase = structure(c(2L, 2L, 3L, 3L, 3L, 1L), .Label = c("G1",
    "G2M", "S"), class = "factor"), old.ident = structure(c(5L,
    5L, 5L, 5L, 5L, 5L), .Label = c("Fibroblast", "T cell", "Macrophage",
    "Stellate", "Acinar", "Endothelial", "Tumor", "B cell", "Mast cell",
    "Ductal", "Islets of Langerhans"), class = "factor")), row.names = c("untreated_TTTGGTTGTCTAATCG-18",
"untreated_TTTGGTTTCCCGAGGT-18", "untreated_TTTGTTGAGAACTGAT-18",
"untreated_TTTGTTGAGCTCGGCT-18", "untreated_TTTGTTGAGTGCCTCG-18",
"untreated_TTTGTTGCACGGTGCT-18"), class = "data.frame")

Been using this code to generate the previous graph.
ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

What is DF ? Can you provide an image of the graph when you are comparing cell cycle in each condition ? Your example "all.combined@metadata" contains only "untreated" conditions, maybe you should update your example to provide few lines of treated conditions too.

geom_bar(stat = "identity") is rarely correct, as bar plots should be used for count-level data only. If you want to use a value, then use geom_col() instead - this is explained in the help for geom_bar.

You don't have to post all the metadata file, but few lines of each condition "treated / untreated" could be helpful.

@dc37 I posted the head and tail of the dataset.

Based on your example, you can use count for each phase of cell cycle for each cell types in each condition. Here, I'm using dplyr package to do that but most likely, you could have the same output using various methods:
DF <- rbind(dfhead, dftail)

library(dplyr)
DF_Count <- DF %>%group_by(orig.ident,Phase,old.ident) %>%
  count() %>%
  ungroup() %>%
  group_by(orig.ident,old.ident) %>%
  mutate(Freq = n/sum(n)*100)

# A tibble: 8 x 5
# Groups:   orig.ident, old.ident [5]
  orig.ident Phase old.ident      n  Freq
  <chr>      <fct> <fct>      <int> <dbl>
1 treated    G1    Fibroblast     1 100
2 treated    G1    Stellate       1 100
3 treated    G2M   Tumor          1  33.3
4 treated    S     Tumor          2  66.7
5 treated    S     Mast cell      1 100
6 untreated  G1    Acinar         1  16.7
7 untreated  G2M   Acinar         2  33.3
8 untreated  S     Acinar         3  50

As you can see, DF_Count has the frequency for each cell types of each phase of the cell cycle in function of the condition. We can use DF_Count to get the following plot. Using facet_wrap, you can create two panels based on the condition column and thus represent the cell cycle of each cell types in function of the treatment condition:
ggplot(DF_Count, aes(x = old.ident, y = Freq, fill = Phase))+
  geom_col()+
  geom_text(aes(label = paste(round(Freq, 2),"%")),position = position_stack(vjust = 0.5))+
  facet_wrap(~orig.ident)

Does it look what you are expecting ?
NB: DF here is the addition of your head and tail of your all.combined dataframe.

How about the other cell types? Can I just make DF <- all.combined@metadata? That way all the cell types are included?

It should work. Here the graph represent only the portion of data you have provided. If you think this answer was useful, do not hesitate to vote or validate it. Learn more here: https://bioinformatics.stackexchange.com/help/someone-answers


text:
stacked bargraph with error bars
I'm attempting to plot a stacked barplot with ggplot2 with this code:
ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

Question: How do I go about adding error bars to the graph?
I tried this code:
ggplot(DF, aes(x = Condition, y = Percent, fill = Cell_Type))+
     geom_bar(stat = "identity")+
     geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5) + geom_errorbar(aes(ymin=Condition, ymax=Condition), width = .2))

But got this error:

Error: Cannot add ggproto objects together.

output of dput(DF)
structure(list(Cell_Type = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), .Label = c("Fibroblast", "T cell",
"Macrophage", "Tumor", "Islets of Langerhans", "Endothelial",
"B cell"), class = "factor"), Condition = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("treated",
"untreated"), class = "factor"), Freq = c(6051L, 1892L, 1133L,
657L, 116L, 868L, 832L, 5331L, 3757L, 1802L, 835L, 287L, 704L,
256L), Percent = c(52.3941466793662, 16.3823707680319, 9.8103731924842,
5.68880422547407, 1.00441596675037, 7.51580223395965, 7.20408693393367,
41.0962072155412, 28.9623805118717, 13.8914585260561, 6.43694110391613,
2.21245760098674, 5.42707369719396, 1.97348134443417)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -14L), groups = structure(list(
    Condition = structure(1:2, .Label = c("treated", "untreated"
    ), class = "factor"), .rows = list(1:7, 8:14)), row.names = c(NA,
-2L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE))

What are len and sd ? can you add a reproducible example of DF, len,sd  ? By curiosity, does your question is related to this one https://bioinformatics.stackexchange.com/questions/11197/stacked-barplot-for-single-cell-analysis?

I didn't change len and sd when I pasted it. I have updated it to Condition and added the output of DF. And yes, it is related to that question. Wasn't sure whether to ask it on the same question.

In your exemple, there is no mention of len and sd, so it is not clear what are these values and how you want to calculate them. Also, does the related question is yours ? because you don't have the same username ?

len and sd are not part of my example. When I pasted the code to the question, I didn't change it. Yes, it is my question, as for the username, I just created a new account as I couldn't log into the other.

Ok, so your question is about to calculate len and sd ? because you don't have these values, right ? For your account, you can see this link: https://bioinformatics.stackexchange.com/help/merging-accounts. Also, when you are registered on one site of StackExchange, you can directly create a new account on another StackExchange website by clicking on topright join Community button, it will automatically link your account on all sites you are registered of.

You are getting Error: Cannot add ggproto objects together. because you lack a ) and geom_text() is not closed. And once this is fixed, you would get Error: Discrete value supplied to continuous scale because geom_errorbar() is for intervals but your Condition variable is a factor.

@dc37 No, my question is how to do a stacked bar graph with the error bars.

@haci How do I change it to an interval?

@user6810 Well, you have to have some min and max values to plot. Depends on why you need these error bars.

I want to indicate the error or uncertainty of each cell type.

If you want to indicate the error or the uncertainty for each cell type, you need to replicate your experiments, ccount each cell type in each replicate, calculate the average percentage of each cell type in each condition and the associated standard deviation.
To answer your question about how to plot a stacked bargraph with error bar, first this option is not integrated into ggplot2 (see this discussion from the creator of ggplot2: https://github.com/tidyverse/ggplot2/issues/1079)
So, you need to calculate the position of each error bar by hand. You can do that by using dplyr. Here, I set the same error bar for everyone as you don't know what is the standard deviation in your data but at least you can observe the structure of the dataframe:
library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  mutate(Cell_Type = factor(Cell_Type, levels = c("B cell","Endothelial","Islets of Langerhans","Tumor","Macrophage","T cell","Fibroblast"))) %>%
  group_by(Condition) %>%
  mutate(SDPos = cumsum(Percent)) %>%
  ggplot(aes(x = Condition, y = Percent, fill = Cell_Type))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))+
  geom_errorbar(aes(ymin = SDPos-SD, ymax = SDPos+SD), width = 0.3, position = "identity")

As you can see, it's pretty ugly... so, if you need to add error bar, I will rather advised you to do a dodged barchart as the following:
library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  ggplot(aes(x = Cell_Type, y = Percent, fill = Condition))+
  geom_bar(stat = "identity", position = position_dodge())+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_dodge(0.9), vjust =-2, size = 3)+
  geom_errorbar(aes(ymin = Percent-SD, ymax = Percent+SD), width = 0.3, position = position_dodge(0.9))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_y_continuous(limits = c(-2,60))

To my opinion, it makes easier the comparison of the distribution of each cell type between conditions.
Hope that it answers your question.

For this scenario, I would advise against a barchart and would use a violin plot or a boxplot instead. The latter two would enable to assess the number of replicates within a group as well as the spread of the replicates.

I agree if you are able to have more than 3 replicates. Otherwise, the spread observed with boxplot or violin plot won't mean a lot. But as it seems that individual values are the count of different cell types from single cell RNA seq, I doubt you can have a lot of replicates. That's why I advised for barchart. But it is also possible to plot individual values and the mean using geom_pointrange.

These numbers and graphs will not mean anything anyway if one does not have many replicates, especially if the data are from human samples and not from a model organism with a rather "clean" genetic background. Moreover, it is even more important to use box/violin plots when the number of replicates are small so that the reader is not mislead about the lack of power in the analyses. I guess, at least hope, that we will not see these plots in the future anyway: https://twitter.com/rafalab/status/1098973538260840449?lang=en

Totally agree ;) and personally, I won't publish those kind of results (at least in this format). But it's well beyond the scope of the question :D

These samples are derived from human samples, and I have replicates for each sample. A total of 32 samples, 16 from each condition. What do you all propose I should use to plot? I'm new at bioinformatics, and assumed a bargraph would be ideal. From reading the comments and answer, I don't have SD for each. I will add that.

If you have 16 for each conditions, you should associate them in a single dataframe and plot a boxplot or a violin plot as advised by @haci.


text:
Filtering paired-end reads with sambamba: avoid discarding reads on the minus strand
I have a BAM file (DNA, shallow whole genome sequencing at ~1X) where I want to filter reads (using sambamba) to keep only those which have a template length > 20 and mapping quality > 20, discarding unmapped reads.
The filter I have used is this:
sambamba view -h \
   -t 5 \
   -F "paired and template_length > 0 and mapping_quality > 20 and not (unmapped or mate_is_unmapped)" \
   -f bam 142-ready.bam \
   -o 142-filter.bam

However, inspection of the reads in IGV led me to think that perhaps reads on the minus strand are being somehow filtered out.
Pre filter:

Post filter:

The only way I don't get such an aggressive filtering is by setting only not (unmapped or mate_is_unmapped). But that, obviously, defeats the purpose as it doesn't include what I need.
I can't, however, see what exactly is wrong with the filter. What is the correct way to filter for template length while not accidentally discarding reads on the minus strand?

The reason is that sambamba does not use the absolute value of the template length, so, if the read maps so the minus strand, the template length is actually negative. Setting a template length > 0 means that all reads mapping to the minus strand are removed.
In fact, you can't still filter for negative template lengths. See this GitHub PR.


text:
umap highlighting two different models
I'm trying to create a umap for single cell data from human samples and ptx samples. I can get the umap to where it shows the umap with the different clusters but I want to show where the ptx samples are and where the human samples are.
My code is as follows:
#create the Seurat object

OD_10K_HUMAN <- CreateSeuratObject(counts = HUMAN_OD_10K.data, min.cells = 0, project = "human")
SD_5K_HUMAN <- CreateSeuratObject(counts = HUMAN_SD_5K.data, min.cells = 0, project = "human")
BNL.5K <- CreateSeuratObject(counts = SD_BNL_5K.data, min.cells = 0, project = "ptx")
BNL.10K <- CreateSeuratObject(counts = OD_BNL_10K.data, min.cells = 0, project = "ptx")
BNM.10K <- CreateSeuratObject(counts = OD_BNM_10K.data, min.cells = 0, project = "ptx")
BNM.5K <- CreateSeuratObject(counts = SD_BNM_5K.data, min.cells = 0, project = "ptx")

#merge data
scData <- merge(BNL.10K, y = c(BNL.5K, BNM.10K, BNM.5K, SD_5K_HUMAN, OD_10K_HUMAN), add.cell.ids = c("A", "B", "C", "D", "E", "F"), project = "HTB2876")

mark the mito genes
mito.genes <- grep(pattern = "^MT-", x = rownames(x = scData), value = TRUE)
length(mito.genes)
scData[["percent.mt"]] <- PercentageFeatureSet(scData, pattern = "^MT-")

scData[["log_nCount_RNA"]] <- log2(scData[["nCount_RNA"]]+1)

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt < 30)
scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))
#DimPlot(scData, reduction = "pca")

numPC = 20
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.4)
scData <- RunUMAP(scData, dims = 1:numPC)

DimPlot(scData, reduction = "umap", label = TRUE)

By default, DimPlot colors the cells by their "identity class", stored in Ident(scData) in your case. This indeed contains the clustering information.
However, you can you use the group.by argument for DimPlot to pass another column from scData@meta.data to color cells by. In your case you would need a metadata column that stores the ptx/human information.
The merge function automatically stores the original "project" information of the individual Seurat objects that you combine under orig.ident in the metadata. Since you set the project names to either "ptx" or "human", the orig.ident column should contain the information you need.
Something like this should work:
DimPlot(scData, reduction = "umap", group.by = "orig.ident", label = TRUE)


text:
Bar plot with error bars
I am trying to plot a box plot/stacked bar graph with error bars.
output of dput(DF)
`structure(list(Cell_Type = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), .Label = c("Fibroblast", "T cell",
"Macrophage", "Tumor", "Islets of Langerhans", "Endothelial",
"B cell"), class = "factor"), Condition = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("treated",
"untreated"), class = "factor"), Freq = c(6051L, 1892L, 1133L,
657L, 116L, 868L, 832L, 5331L, 3757L, 1802L, 835L, 287L, 704L,
256L), Percent = c(52.3941466793662, 16.3823707680319, 9.8103731924842,
5.68880422547407, 1.00441596675037, 7.51580223395965, 7.20408693393367,
41.0962072155412, 28.9623805118717, 13.8914585260561, 6.43694110391613,
2.21245760098674, 5.42707369719396, 1.97348134443417)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -14L), groups = structure(list(
    Condition = structure(1:2, .Label = c("treated", "untreated"
    ), class = "factor"), .rows = list(1:7, 8:14)), row.names = c(NA,
-2L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE))`

Gives me the following dataset.
       Var1      Var2 Freq
1   Fibroblast   treated 6051
2       T cell   treated 1892
3   Macrophage   treated 1133
4     Stellate   treated  893
5       Acinar   treated  148
6  Endothelial   treated  868
7   Fibroblast untreated 5331
8       T cell untreated 3757
9   Macrophage untreated 1802
10    Stellate untreated 1061
11      Acinar untreated 2786
12 Endothelial untreated  704

library(dplyr)
DF <- df %>% rename(Cell_Type = Var1, Condition = Var2) %>%
  group_by(Condition) %>%
  mutate(Percent = Freq / sum(Freq)*100)

# A tibble: 12 x 4
# Groups:   Condition [2]
   Cell_Type   Condition  Freq Percent
   <fct>       <fct>     <int>   <dbl>
 1 Fibroblast  treated    6051   55.1
 2 T cell      treated    1892   17.2
 3 Macrophage  treated    1133   10.3
 4 Stellate    treated     893    8.13
 5 Acinar      treated     148    1.35
 6 Endothelial treated     868    7.90
 7 Fibroblast  untreated  5331   34.5
 8 T cell      untreated  3757   24.3
 9 Macrophage  untreated  1802   11.7
10 Stellate    untreated  1061    6.87
11 Acinar      untreated  2786   18.0
12 Endothelial untreated   704    4.56

Calculated SD/Mean with this:
my_sum_1 <- DF %>%
    group_by(Cell_Type) %>%
    summarise(
        n=n(),
        mean=mean(Freq),
        sd=sd(Freq)
    ) %>%
    mutate( se=sd/sqrt(n))  %>%
    mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

Then plotted with this:
ggplot(my_sum_1) +
    geom_bar( aes(x=Cell_Type, y=mean), stat="identity", fill="forestgreen", alpha=0.5) +
    geom_errorbar( aes(x=Cell_Type, ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
    ggtitle("Cell Types")

I want to plot based on condition. Not sure how I would do that. I tried the fill command but doesn't seem to work as the object is not found, I can't seem to add it back to the data set. How do I include the condition in the bar chart?
When I'm calculating sd, I think it's combining the cell type value into one, rather than taking it individually. How do I change that?
dput(head(all@meta.data))
structure(list(orig.ident = c("PDAC", "PDAC", "PDAC", "PDAC",
"PDAC", "PDAC"), nCount_RNA = c(7945, 7616, 7849, 5499, 853,
1039), nFeature_RNA = c(2497L, 2272L, 2303L, 2229L, 509L, 588L
), percent.mt = c(2.63058527375708, 3.7421218487395, 4.9433048796025,
0.490998363338789, 0.234466588511137, 0.192492781520693), RNA_snn_res.0.5 = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), seurat_clusters = structure(c(5L,
5L, 5L, 5L, 1L, 6L), .Label = c("0", "1", "2", "3", "4", "5",
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16",
"17", "18", "19"), class = "factor"), S.Score = c(0.0171539298241524,
-0.00599651468836135, 0.011058574403656, -0.0424413740104151,
0.03549521163095, -0.0243347616753686), G2M.Score = c(-0.0591387992008088,
0.0253275642795205, -0.0512402869839816, -0.0239076248512967,
-0.0263164126515597, -0.0339086517371782), Phase = structure(c(3L,
2L, 3L, 1L, 3L, 1L), .Label = c("G1", "G2M", "S"), class = "factor"),
    old.ident = structure(c(1L, 1L, 1L, 1L, 1L, 3L), .Label = c("Fibroblast",
    "T Cell", "Endothelial", "Tumor", "Stellate", "Macrophage",
    "B Cell", "Mast Cell", "Acinar", "Endocrine", "Exocrine"), class = "factor")), row.names = c("PDAC_AAACCCAGTCGGCTAC-1",
"PDAC_AAACGAAGTCCAGGTC-1", "PDAC_AAAGAACAGCAAGTGC-1", "PDAC_AAAGAACAGGATGAGA-1",
"PDAC_AAAGTGACATCCGTTC-1", "PDAC_AACAACCAGGAAGTAG-1"), class = "data.frame")

Your issue is coming by the fact that when you calculate mean, se, ... you removed the Condition column, here is your my_sum_1 results:
# A tibble: 7 x 6
  Cell_Type                n  mean    sd    se     ic
  <fct>                <int> <dbl> <dbl> <dbl>  <dbl>
1 Fibroblast               2 5691   509. 360    4574.
2 T cell                   2 2824. 1319. 932.  11849.
3 Macrophage               2 1468.  473. 334.   4250.
4 Tumor                    2  746   126.  89.0  1131.
5 Islets of Langerhans     2  202.  121.  85.5  1086.
6 Endothelial              2  786   116.  82    1042.
7 B cell                   2  544   407. 288    3659.

To calculate the mean, sd, ... from each cell types in each conditions, you need to group by Cell_Type and Condition by writing:
my_sum_1 <- DF %>%
  group_by(Cell_Type, Condition) %>%
  summarise(
    n=n(),
    mean=mean(Freq),
    sd=sd(Freq)
  ) %>%
  mutate( se=sd/sqrt(n))  %>%
  mutate( ic=se * qt((1-0.05)/2 + .5, n-1))

But now, your table looks like:
# A tibble: 14 x 7
# Groups:   Cell_Type [7]
   Cell_Type            Condition     n  mean    sd    se    ic
   <fct>                <fct>     <int> <dbl> <dbl> <dbl> <dbl>
 1 Fibroblast           treated       1  6051    NA    NA    NA
 2 Fibroblast           untreated     1  5331    NA    NA    NA
 3 T cell               treated       1  1892    NA    NA    NA
 4 T cell               untreated     1  3757    NA    NA    NA
 5 Macrophage           treated       1  1133    NA    NA    NA
 6 Macrophage           untreated     1  1802    NA    NA    NA
 7 Tumor                treated       1   657    NA    NA    NA
 8 Tumor                untreated     1   835    NA    NA    NA
 9 Islets of Langerhans treated       1   116    NA    NA    NA
10 Islets of Langerhans untreated     1   287    NA    NA    NA
11 Endothelial          treated       1   868    NA    NA    NA
12 Endothelial          untreated     1   704    NA    NA    NA
13 B cell               treated       1   832    NA    NA    NA
14 B cell               untreated     1   256    NA    NA    NA

Because in your initial dataset, you only have one value for each cell type in each condition, so you can calculate sd, sem. You need to first merge your replicates into your dataframe in order to be able to calculate your sem.
Then, you will be able to do:
library(ggplot2)
ggplot(my_sum_1, aes(x=Cell_Type, y=mean, fill = Condition)) +
  geom_col(alpha=0.5) +
  geom_errorbar( aes(ymin=mean-ic, ymax=mean+ic), width=0.4, colour="orange", alpha=0.9, size=1.5) +
  ggtitle("Cell Types")

Hope it helps you to figure it out what is your issue with your current data.

How do I merge a dataframe? These aren't replicates though. For each cell type, I have two conditions.

If you made replicates of your experiment, you should have multiple tables similar to  DF (one for each replicate) with the count of each cell types per conditions. If so, you can merge all of them (take a look for function ?merge() or ?full_join() from dplyr (https://dplyr.tidyverse.org/reference/join.html)). If all of your tables have the same format and colnames, you can consider using rbind (and adding first a condition column to each table corresponding to the replicate number).

I uploaded the question with the dataset before I calculated the frequencies and percentages.

I don't think it changes anything. This table does not represent replicates right ? Do you have replicates for this experiment ? What I mean by replicates is that you run a scRNAseq several times on your biological samples, and you identify the percentage of cell types in each of this runs. How do you design your experiments ? Can you describe it a little bit ? (no need to know exact parameters, just what you are trying to do and what materials did you use at the beginning)

I run one sample on the 10x chromium box. This data set contains samples from treated and untreated patients. 15 treated and 16 untreated samples. I run all the raw files through Seurat.

I uploaded the metadata from Seurat.


text:
What genes are driving a certain pathway
I've currently run pathway analysis on a set of genes. I'm now trying to figure out how to see which genes are driving each pathway.
cleaned_ductal_markers <- tumor_markers
markers15  = cleaned_ductal_markers %>% filter(cluster=="PSEUDOPAPILLARY", avg_logFC>0.5, p_val_adj<0.01)
markers15

#head(genes_map)
markers_entrez = genes_map[which(genes_map$SYMBOL %in% markers15$gene),"ENTREZID"]
kk <- enrichKEGG(gene         = markers_entrez,
                universe      = genes_map[, "ENTREZID"],
                 organism     = 'hsa',
                 pvalueCutoff = 0.05)
#head(kk)
dotplot(kk, showCategory=10)
options(repr.plot.width=7, repr.plot.height=5)
emapplot(kk, showCategory=12, font.size =16) + set_graph_style(family = "Arial Narrow", face = "plain",  text_size = 15, text_colour = "black")

I'm now given a list of pathways but now I'm trying to see which genes are driving those pathways.
How do I see the genes?

What’s your actual question?

@user438383 My question is how do I see which genes are driving those pathways?


text:
ViolinPlot x-axis
I'm trying to set the x-axis in my Violin Plot to go by each patient rather than it being scattered.
VlnPlot(obj, features = c("classical_basal"), pt.size = 0) + NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

The above code gives me the x-axis but there is no pattern or order.
I would like the x-axis to be grouped with the HTB, for example ALL of the HTB2867 first, followed by HTB2876, etc...

Seurat does violin plots using ggplot, which orders based on the factor levels. You need to re-factor your grouping category so that the levels matches your preferred order:
# default level sort order is alphabetical
obj[["grouping"]] <- factor(unlist(obj[["grouping"]])) # refactor
VlnPlot(obj, group.by = "grouping", features = c("classical_basal"),
        pt.size = 0) +
  NoLegend() + theme(axis.text.x = element_text(angle = 90, hjust = 1))

How do I re-factor the grouping category?

It is difficult to say exactly what you would like to achieve without an example.
But if you would like an individual axis based on a factor like patient maybe you want to use facet_wrap for the Patient variable like in the example below?
library(ggplot2)

data <- iris

ggplot(data, aes(x=Sepal.Width, y=Sepal.Length, fill=Species)) +
  geom_violin() +
  facet_wrap(~Species, scales="free_x")
```

I have attached an image as well as a better description.


text:
Creating new seurat object with new matrix
I'm trying to do a cross-species comparison between the patient TME vs. the PTX TME to understand gene expression conservation. I have already converted between mouse and human genes giving me a resultant gene expression matrix that is just one set of (homologous) genes X cells for all mouse and human cells.
musGenes <- rownames(mouse_only)
# Basic function to convert mouse to human gene names
convertMouseGeneList <- function(x){
require("biomaRt")
h = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
m = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
genesV2 = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = x , mart = m, attributesL = c("hgnc_symbol"), martL = h, uniqueRows=T)
humanx <- unique(genesV2[, 2])
# Print the first 6 genes found to the screen
print(head(humanx))
return(humanx)
}
genes <- convertMouseGeneList(musGenes)

Now how do I create a new Seurat object with the converted mouse gene list, such as rename the matrix, so that the rownames/colnames are human symbols, rather than mouse symbols?

is this thread useful? They do rename rownames of the just created seurat object.


text:
Plotting infercnv results
I'm working with matched single cell data, where we have treated and untreated samples for the same patient. I ran CNV analysis using the infercnv package.
I've followed the tutorial:
# data matrix
counts_matrix <- scData@assays$RNA@counts
meta = data.frame(labels = Idents(scData), row.names = names(Idents(scData)))
unique(meta$labels) # check the cell labels

# Create the infercnv object
infercnv_obj <- CreateInfercnvObject(raw_counts_matrix = counts_matrix,
                                     annotations_file = meta,
                                     delim = "\t",
                                     gene_order_file = "hg38_gencode_v27.txt",
                                     ref_group_names = NULL)

# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir="output_dir",  # dir is auto-created for storing outputs
                             cluster_by_groups=T,   # cluster
                             denoise=T,
                             HMM=T
                             )

Now, I'm not sure how to plot the results. Any help would be appreciated.
When running the below recommendation I get this:
INFO [2023-05-19 12:49:22] ::plot_cnv:Start
INFO [2023-05-19 12:49:22] ::plot_cnv:Current data dimensions (r,c)=34108,17624 Total=82463956 Min=0 Max=55099.
INFO [2023-05-19 12:49:22] ::plot_cnv:Depending on the size of the matrix this may take a moment.
<sparse>[ <logic> ]: .M.sub.i.logical() maybe inefficient
Warning: sparse->dense coercion: allocating vector of size 4.5 GiBINFO [2023-05-19 12:50:12] plot_cnv(): auto thresholding at: (-2.725632 , 3.000000)
INFO [2023-05-19 12:50:35] plot_cnv_observation:Start
INFO [2023-05-19 12:50:35] Observation data size: Cells= 17624 Genes= 34108
INFO [2023-05-19 12:50:35] clustering observations via method: ward.D
INFO [2023-05-19 12:50:35] Number of cells in group(1) is 331
INFO [2023-05-19 12:50:36] group size being clustered:  331,34108
Error in parallelDist(data_to_cluster, threads = infercnv.env$GLOBAL_NUM_THREADS) :
  x must be a matrix or a list of matrices.

If you haven't already, I recommend checking out the infercnv YouTube tutorial (https://github.com/broadinstitute/inferCNV/wiki/Running-InferCNV) for detailed instructions on running InferCNV. Additionally, you can explore Uphyloplot2 (https://github.com/harbourlab/uphyloplot2) for building a tree from the output generated by InferCNV.
I noticed that you didn't use the untreated sample as the normal reference for your CNV analysis. I'm curious to understand your reasoning behind this. In my case, I also have an untreated sample that isn't diploid, along with two treated samples. Initially, I planned to use the untreated sample as the "normal" reference. However, considering that InferCNV assumes a normal diploid reference, I'm now uncertain if using the untreated control would be appropriate.
Please let me know your thoughts, any help would be much appreciated!

Personally I would not use intercnv_obj throughout, because is there was a bug in line 2, i.e. using = instead of <- then the created intercnv_obj might be trying to be passed into plot_cnv.
Could you try
data(infercnv_obj)

plot_cnv(
  infercnv_obj,
  out_dir = ".",
  title = "inferCNV",
  obs_title = "Observations (Cells)",
  ref_title = "References (Cells)",
  cluster_by_groups = TRUE,
  cluster_references = TRUE,
  plot_chr_scale = FALSE,
  chr_lengths = NULL,
  k_obs_groups = 1,
  contig_cex = 1,
  x.center = mean(infercnv_obj@expr.data),
  x.range = "auto",
  hclust_method = "ward.D",
  custom_color_pal = NULL,
  color_safe_pal = FALSE,
  output_filename = "infercnv",
  output_format = "png",
  png_res = 300,
  dynamic_resize = 0,
  ref_contig = NULL,
  write_expr_matrix = FALSE,
  write_phylo = FALSE,
  useRaster = TRUE
)

If the above doesn't work ...
What I think is happening is the input data is the problem. This precise error is already known and was reported on Broad's GitHub here:
https://github.com/broadinstitute/infercnv/issues/415
The responder (Broad Institute) believed it was the version problem, i.e. the version was outdated. However a separate respondee believed it was the input data format and suggested shifting away from a sparse matrix. That would make sense in context to the error being thrown.

Hi , I get an error when running that. I updated my question with the message.

@mmpp I have modified the code and I've included an error log on GitHub. Apologies for the delay. If problem continued I'd recommend referring to this GitHub error log ... https://github.com/broadinstitute/infercnv/issues/415 or else filing a separate error and referring to this error log. It is not clear whether Broad has fully resolved this - but they already know about it. Good luck :-)

What did you change?

data(infercnv_obj) prior plot_cnv

It’s absolutely fine to use = as an assignment operator in R, it’s largely a myth that <- is ‘better’ in some way. It’s certainly not a ‘bug’ by any means.


text:
functions in seurat to calculate the gene count per cluster
I'm running FindAllMarkers
markers <- FindAllMarkers(scData_scTransform_PTX_only_res04_20PC_6000g_30mt_sb_removal_20210713_ann_exp0_counts, min.pct = 0.25, logfc.threshold = 0.25)

But was wondering if there is a function in seurat to get the raw (non-normalized) expression data from each cell type cluster

Hi, could you please add more context to this question so that we can understand why you want to find out this information? I've calculated cell counts per cluster, and visualised gene counts per cluster using scatter plots, but haven't yet run into a case where I'd need to work out gene count per cluster as a single statistic (whatever that means).

@mmpp could it be that you meant to compare expression profiles of some genes (by means of a boxplot, for instance), to evaluate the differences between two conditions, in the same cell type, i.e. a cluster?

are you trying to calculate pseudobulk average?