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

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

Yes, this is expected. I am not sure what you mean by "reference" vs "database", but what you do is to index a fasta file or a collection of fasta files. That index you create by that is what bowtie2 alignes against. It does not know whether the index consisted of one, ten or a thousand fasta files, it just maps against the entire index and then check whether it can uniquely (or not) place the reads. If you combined two genomes into one index and the genomes are very similar then indeed many reads will map equally well to them and that manifests as MAPQ 0 (or low).

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


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

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

source: Stack Exchange

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

sort: votes


text:
Error in DESeq2 analysis with dimnames
I'm getting the following error in a DESeq2 analysis and I can't seem to figure out the issue:
Error in `assays<-`(`*tmp*`, withDimnames = withDimnames, ..., value = `*vtmp*`) :  please use 'assay(x, withDimnames=FALSE)) <- value' or 'assays(x, withDimnames=FALSE)) <- value' when the dimnames on the supplied assay(s) are not identical to the dimnames on RangedSummarizedExperiment object 'x'
Here is my code:
library(DESeq2)
library(edgeR)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))

samplenames <- substring(colnames(x), 12, nchar(colnames(x)))

colnames(x) <- samplenames
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
                     "Basal", "ML", "LP"))
x$samples$group <- group

#DESeq2 Analysis
phenodata <- new("AnnotatedDataFrame", data = data.frame(celltype = x$samples$group))
counts <- x$counts
ids <- c(1,3,4,6,7,9)

tmpcounts = counts[, ids]
tmpph = as(phenodata[ids, ], "data.frame")
rownames(tmpph) <- colnames(tmpcounts)

dds = DESeqDataSetFromMatrix(countData = tmpcounts, colData = tmpph, design = ~ tmpph$celltype)

Here are my package versions
> packageVersion("DESeq2")
[1] ‘1.30.1’

> packageVersion("SummarizedExperiment")
[1] ‘1.20.0’

What works is:
dds = DESeqDataSetFromMatrix(countData = data.frame(tmpcounts), colData = tmpph, design = ~celltype)

Mind the data.frame() around the count data object (which was a matrix) and removal of tmpph$ for the design because design should reference one or more columns of the colData. In newer versions of DESeq2 that by the way works without the data.frame() but I am not sure why. You might want to dig that out yourself, looking at the source code of how exactly the package (SummarizedExperiment) feeds data as an assay into a SummarizedExperiments object.

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


text:
scRNA-Seq pseudobulk for marker detection
A standard approach for scRNA-Seq is to partition the single cells into individual clusters, then use a Wilcoxon test to find markers that characterize each cluster (or other statistical methods that consider single cells as replicates).
Recently, there has been a push towards using pseudobulk approaches for "differential state" (DS) analysis, which avoid the pseudoreplication bias of considering single cells as replicates. See e.g. Crowell (2020), Zimmerman (2021), Thurman (2021), Squair (2021), or the OSCA book.
However, all the above papers seem to use pseudobulk methods for DS of a single cluster between two conditions (e.g. drug treatment), or occasionally to compare a specific pair of clusters. I haven't seen any example of using pseudobulk to find markers distinguishing clusters (comparisons of one cluster against all others).
I can think of a few reasons, for example:

these pseudobulk methods can be more computationally intensive,
perhaps the Wilcoxon test is good enough to simply find markers,
pseudobulk methods can only work when each cluster is represented in several samples, which is hard to guarantee for every single cluster in the dataset.

On the other hand it seems to me a pseudobulk approach should yield more relevant markers (genes that are robustly enriched across biological replicates). So, am I missing something? Is there a statistical reason not to use pseudobulk approaches for finding cluster-specific markers?

The only reason I see for not using pseudobulks is lack of biological replication. Even though single-cell sequencing is common these days and standardized it is still expensive. If you have several conditions plus biological replicates and need many cells, maybe due to heterogeneity, the bill can easily reach the 5-digit realm of money. Biological replication is always desired but sometimes simply not possible. We recently did a study in which we were interested in organ heterogeneity in different conditions and both mouse work and library prep was so extensive that even hash tagging the mice we had to pool was simply not possible. So we were bound to using each cell as a replicate. I prefer limma-trend for this, not the Wilcox test. It is not computationally very demanding, neither is pseudobulks. But even if it was, I would happily wait 30min rather than 10sec to get the best possible rather than the fastest results. Pseudobulks conceptually have the advantage that you basically are forced to have biological replication (because only then they make sense) so findings are intrinsically more robust. Theoretically, you can of course have biological replicates and still treat each cell as replicate, and then blocking or adding the biological replication information into the DE testing design, but the statistical rigor of methods used for pseudobulks (DESeq2, edgeR and other) is for me the reason to always use pseudobulks if I can.
I personally use pseudobulks also for markers. After all marker detection comes downstream of DE testing by aggregating the results into lists that contain unique or enriched genes per cluster. I see no reason to use multiple testing strategies.

If you want a statistical reason not to use pseudobulk to do DE for marker finding, I think it is the same reason not to use any DE method for marker finding: the double dipping problem (see Gao and Witten 2020).
Briefly, if you cluster a dataset based on genes (we can safely say that clustering in PCA space is still based on genes), then asking which genes are differentially expressed between clusters is a statistical fallacy.
We could just as easily divide apples in 2 groups according to their shade of red and then ask which RGB channels exhibit a statistically significant difference.
This means that even in the absence of real structure, you will identify very low p-values for your DE test. The OSCA book section you mention links to another section in the advanced book that tells you the same thing with a small practical demonstration generating a 200-dimensional Gaussian and clustering it, then testing for differences in clusters.
This does not mean, however, that markers found this way are meaningless; it simply means that you cannot use a p-value cutoff to define whether a gene is a marker in this scenario, or say that since you have statistically significant markers two clusters are different.
You can still rank these markers by fold change, difference in percentage expression, AUC, entropy, etc as well as p-value. Moreover, markers identified using the traditional methods tend to make biological sense, and statistically rigorous methods are still in development (there was a beautiful talk by Anna Neufeld on an approach involving count splitting, you can read the preprint here).
The test you suggest is interesting, but I don't think you would solve the double dipping problem; this being said it looks more robust and appropriate than using every single cell as a replicate.

Thank you for the useful addition, indeed the pseudobulk approach does not solve the double-dipping problem, and will suffer from the same p-value inflation as a Wilcoxon on the single cells.


text:
Select synonymous sites from a multiple sequence alignment
Could someone kindly recommend a tool or R package that can identify synonymous sites in a multiple sequence alignment?
I wish to select those taxa for tree reconstruction and other downstream analyses within subset of the alignment.
To respond to the questions below: yes, I am interested in identifying synonymous sites within a multiple sequence alignment (DNA) of ~70 virus samples with limited genetic diversity. My goal is to subset the multiple sequence alignment including synonymous sites only, because we believe that there are parallel non-synonymous mutations occurring that do not provide information on shared evolutionary history.
Thank you @M_ for the really helpful information below. In Mega, the S tool identifies singleton sites rather than synonymous sites in a multiple sequence alignment. I'm wondering if there is another tool or a PAML option to subset the synonymous sites?
Thank you in advance and best!

Do we also need to handle the translation, or is it already a protein alignment? Are we talking about two sequences or more? What kind of output do you need? Visual? Or text? Something else?

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

I proposed Syn-SCAN (below). This was used for HIV and as you are looking at viruses seems a good fit.
There's a simpler way ... just make a invariable amino acid alignment, i.e. discard all amino acid variation. MEGA will do that, or use an algorithm. You then simply check the remaining nucleotide sequences for identity and discard the identity (lots of algorithms do this).
What you are left with is, is what you seek, PURE SYNONYMOUS mutations.

discard taxa with amino acid variation
retain taxa nucleotides with variation.

However, just discarding amino acid variation and retaining everything else will either by identity or synonymous. To be honest, its a virus so the amount of sequence identity between taxa is likely to be very little for protein loci.
Its really that simple. I honestly wouldn't over complicate it.
Sorry for the singleton vs. synonymous error, I haven't looked at MEGA in ages. Its a bit bizarre it discarded the only thing it was actually used for.

The short answer to the question is the 3rd codon position is the traditional approximation for synonymous mutation. A small proportion of 1st codon position sites are neutral and a small proportion of 3rd codon position are non-synonymous.
I agree with @terdon about and @Community regarding additional information sought, because synonymous mutations are heavily subject to compositional bias and that is a big subject. The key issue is the whether the taxa/species in the alignments are known for compositional bias and the genetic distance between different members of the alignment.
The traditional package used for synonymous mutations is/was MEGA, currently version 11. Albeit, recent 'upgrades' seem to remove the synonymous distance matrix calculation. HOWEVER, they may have transferred this capability onto its what appears to be upgraded GUI alignment viewer. What I think is if you:

load the alignment into the GUI;
click alignment (big button far left of menu);
you will now see the alignment;
in the middle of the menu bar are 4 buttons 'C', 'V', 'Pi' and 'S'
I think 'S' = synonymous
All synonymous sites in the alignment are now highlighted and there is likely a method of extracting this information.

C=constant, V=variable Pi=Nei approximation of diversity (population genetics)
If this is correct then those are the subset sought, please note however you will need to select/deselect the taxa manually and reassess the synonymous sites after each 'point and click' iteration (the S button will remain active I suspect between each selection/deselection). What you appear to want is a synonymous tree and this function I believe was removed from MEGA.
MEGA is rarely used except for teaching because its functionality is available elsewhere that permit a given algorithm to be integrated into a data pipeline, which MEGA would not permit (might have recently changed). Hence my recollection of MEGA is a little vague.

The package used for synonymous work is PAML and remains the undisputed standard without question in the context of formal selection analysis, which is heavily dependent on a estimate of synonymous mutations. What you refer to as 'synonymous mutation' isn't the same as the phylogenetic definition, i.e. your definition depends on the other members of the alignment, whereas this isn't an issue in the phylogenetic definition (because the tree circumvents this, i.e. delineates subsets) except if compositional bias occurs ... huge amounts of work is/has been dedicated to compositional bias.
Final point is the widely used phylogenetic library in R called ape will not perform a synonymous mutation calculation - I think, however phangorn library in R might do, but is more likely within the context of selection analysis. Thus to answer the question directly, I don't think R will do this calculation.

It transpires S= singleton (sorry) ....
Have a look at Syn-SCAN. The link is here https://hivdb.stanford.edu/pages/synscan.html . I have never used it and the download link is on the website. The good thing is its written in Perl so it would be easy to modify and check (if you can read/write Perl).
PAML: you can't detach the synonymous mutation calculation from the non-synonymous mutation calculation, without delving into the C-code (which I wouldn't recommend).

"Final point is whether R's ape": ape? Should that be API or something else entirely?

Thank you @M_ for the above! I edited my question in response. Thanks again!

Hi @terdon .. sorry ape is a popular phylogenetics library in R. Personally, I think its over-used. Its not an API ... I will clarify the text.

Thanks! As for "Personally, I think its over-used", that's how I feel about R :P


text:
Imputing small region of the genome
If I'm looking for a specific SNP in my SNP-Chip data and it isn't there, are there any tools that let me quickly impute that SNP from surrounding SNPs rather than running a lengthy 'whole chromosome' imputation job?
If so, roughly how many upstream and downstream SNPs are generally required to predict the haploblock 'accurately'?

I would take something like 1.5Mb either side of the SNP of interest (so a 3Mb chunk), subset that out from the target and reference panel and then perform imputation on that chunk alone. I think that size will give you more than enough LD, since LD decays by about 500kB in most human populations by then.
shapeit4 has an option to just impute a single region:
shapeit4 \
    --input unphased.bcf \
    --map chr20.b37.gmap.gz \
    --region chr21:1231213-1413143 \ ## change this bit
    --reference reference.bcf \
    --output phased.bcf

Thanks, that's clear enough. I'm just looking for a tool to make that easy for me ;-)

@dan see edit, shapeit4 can impute just a particular region


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

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

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

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

Adding my answer from the crosspost:
I would not use the baseMean for any filtering as it is (at least to me) hard to deconvolute. You do not know why the baseMean is low, either because there is no difference between groups and the gene is just lowly-expressed (and/or short), or it is moderately expressed in one but off in the other. The baseMean could be the same in these two scenarios. If you filter I would do it on the counts. So you could say that all or a fraction of samples of at least one group must have 10 or more counts. That will ensure that you remove genes that have many low counts or zeros across the groups rather than nested by group, the latter would be a good DE candidate so it should not be removed. Or you do that automated, e.g. using the edgeR function filterByExpr.


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

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

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

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

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

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

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

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


text:
Pymol: select low confidence regions from AlphaFold pdb file
I have downloaded a predicted structure from AlphaFold as a pdb file (https://alphafold.com/entry/O75376) and loaded it into Pymol (2.3.0). There is quite a large portion of the structure that was modelled with very low confidence (pLDDT < 50). Is there a way to select these regions in the structure so that I can assign a different representation to them? Thank you.

Answer from @matteo-ferla, converted from comment:
They are stored as b-factors. This is well described in the PyMol wiki.
[please edit to improve this answer, if possible]

use:
Spectrum b, Red_white_blue

It is not clear what this means.


text:
Repeating indexes of atoms in PDB file
I have a PDB file in which the indexes of atoms are repeated (or better say refreshed) every 99999.
I am using MDtraj software to analyze this pdf file, so these repeating indexes cause some problems.
Does everyone encounter this problem and how can I fix it?

Inconsistencies like randomly repeated atom indices, maxed out indices at 99999, integer overflowing indices to the negatives (yup, there's a program that has this bug), bad spacing, jumbled order etc. in a PDB are tolerated differently in different programs. Some will carp at missing atoms or non-standard atom names, which are not formatting errors.
In terms of atom indices, these are used in CONECT entries, but several programs bond by proximity or by known topology —CONECT entries are absent for canonical amino acids. Repeated residue indices and chain letter repeats cause issues with LINK and often with peptide chain sequence.
However, in a multimodel PDB the numbers in each model are expected to be reset at every MODEL so if the atom indices reset after the end of a model this is 100% expected.


text:
What is the difference between Normalized Expression in EdgeR vs DESeq2?
I am trying to access the normalized expression in both edgeR and DESeq2, yet the results are different. Does anyone know why?
How to get normalized expression using edgeR:

y <- DGEList(...)
y <- calcNormFactors(y, method = "TMM")
v <- voom(y, design, plot=FALSE)

# normalized expression
head(v$E)

How to get normalized expression using DESeq2:

dds <- DESeqDataSetFromMatrix(...)
dds <- estimateSizeFactors(dds)

# normalized expression
normalized_counts <- counts(dds, normalized = TRUE)
head(normalized_counts)

Anyone know why the results are different? Aren't they both doing a version of log(cpm) under the hood? Or am I not comparing the right results from DESeq2?

First, and most obviously, two independent methods will never fully agree.
Second, while the DESeq2-related code in fact gets the normalized counts based on the size factors the edgeR/voom code extract the voom-transformed counts, but that is not the normalized counts itself. I would need to go back to the limma docs and paper to check what exactly the E slot stores, but I think (please edit for a definite answer) it's normalized counts after the voom transformation. In contrast, the edgeR normalized counts in the sense of what woul be most comparable to the DESeq2 normalized counts would be obtained with the edgeR::cpm() function. Check the manual. At no point of either the limma or edgeR manual it says that v$E should be used for anything outside of the linear modelling. It clearly recommends to use cpm(). Limma and edgeR are packages that have hundreds of forum questions and blog entries, many functions do exist and the packages change over time, so indeed it can be confusing what to do, but therefore I strongly recommend to consult the current manual at Bioconductor as the sole reliable reference on how to get values from the objects and how to run your analysis.
By the way, don't be surprised that the DESeq2 and edgeR normalized counts will be on different scales. DESeq2 gets the counts by dividing the raw counts by the size factors while edgeR first calculates naive per-million counts and then corrects these with the TMM factors calculated by calcNormFactors, so the magnitude of the counts can be quite different, simply by how it is calculated.


text:
Extract chromosome arms' coordinates for specific reference
I would like to get chromosome arms' coordinates for genome reference hs37d5 (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/), but I could not find them either online, or the directory listed above.
Is there a way to extract such coordinates from one of the available files (e.g., from the fasta)?

The chromosome sizes can be retrieved with (to use for sanity check):
> ./faSize -veryDetailed -tab hs37d5.fa.gz
1   249250621   23970000    225280621   225280621   0
2   243199373   4994855 238204518   238204518   0
3   198022430   3225295 194797135   194797135   0
4   191154276   3492600 187661676   187661676   0
...

(first column is chromosome ID, second column is chromosome size).
Afterwards, the centromere coordinates are retrieved through (in this case I used hg19 since I'm working with hs37d5):
> curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz" \
  | gunzip -c | grep acen | sed 's/chr//g' | sort -n > hg19.centromeres

The coordinates for each chromosome arm can be computed from the chromosome size and centromere coordinates.
N.B. As per UCSC documentation, cytoband's start is 0-based, while end is 1-based.
As a sanity check you might want to compare the end of the last cytoband for each chromosome, with the chromosome size (they should be the same, since end is 1-based).


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:
Krakenuniq mmap error
I built Krakenuniq on an HPC, and I've been trying to run it on my data using the following command:
krakenuniq --db DBDIR --threads 1 All_reads.fastq

This definitely isn't powerful enough to get any reasonable output, but this is just a test run. However, I got the following error:
classify: unable to mmap ./DBDIR/database.kdb: Invalid argument

I haven't been able to find anything about this error online. DBDIR is 1TB, which seemed huge, so I was wondering if there might have been a problem with the build, but I saw no issues in the log.

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

Is DBDIR in your working directory? Did you perhaps mean to write --db ${DBDIR}?

That's not how it's done in the documentation, and in any case it causes its own error:

krakenuniq: unable to find --threads in $KRAKEN_DB_PATH (undefined) DBDIR appears to be in the working directory (ls lists it there).

I have run into the same issue with the default Kraken using a prebuilt database at the university HPC.

@schmiggle Did you ever find a solution to this? I can see the database file (for example, ./DBDIR/database.kdb) is there, but I get the same error.

I switched to using Kraken2 instead because of issues like this. Didn't ever find a solution.


text:
Phylogenetic tree rooting in shotgun metagenomics
But I have some weeks fighting with this issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account tree distance branches metrics.
I have one tree for Archea and another tree for bacteria (both unrooted), created from IQTREE from the protein alignments obtained by GTDB-Tk software. I tried to use the software Root-Digger, but only works with nucleotide alignments and not protein alignments. Do you know another software capable (working with protein alignments) to infer more objectively the root rather than define arbitrarily a root through ape R package?

Root-digger is described as "a non-reversible Markov model", which would be assumed to be a directional mutation matrix, thus T->A would have a different probability to A->T using Bayesian approach to branch-length estimation
IQ-TREE is an ultra-rapid method for non-parametric bootstrapping which can be essential for extremely large phylogenies.

Hi @MagiBC this really two questions. Could you kindly focus on your main question and submit the second question as a separate request?

Hi @M__, yes sorry are two questions. Ok, thank you and I take into account your suggestion. The second question is at: https://bioinformatics.stackexchange.com/questions/19018/phylogenetic-tree-building-from-progenomes-database-for-shotgun-metagenomics

One option is to use rootstrap from IQTREE version newer than 2.1.3 (http://www.iqtree.org/doc/Rootstrap). However for at the moment I am not capable to run the first toy tutorial example in the abovementioned question. It's not a stable version. I write in google groups of IQTREE. I keep you informed if I had news.

The answer really is that a root is determined using external evidence, rather than from the internal branches.
An algorithm such as root digger can be used to estimate the root placement based on a branch that for example might minimise the anomalous branch lengths in all other taxa. Excessive comparative branch-lengths could indicate an error in rooting the tree and this can be the basis of the estimate. However, it is an estimate.
R libraries such as ape are more likely to simply mid-point root the tree, which isn't really cool, but its easy to calculate.
The three approaches for rooting are either:

Using a known outgroup to the taxa of interest, which unequivocally predates their divergence.
External information such the evolution of key morphological features of the organism.
Published studies.

Point 1 This can be used singly to find the root - a "root-finding tree" - and the tree of interest presented by constructing a second tree without this root in the alignment. For example the Archaea predate bacteria, one tree determines where the Archaea attach to taxa of interest. If the outgroup is very distinct this could cause artefacts.
Point 2 Finding an outgroup much close to the taxa of interest using external information, for example gram-negative versus gram-positive bacteria, e.g. if all the bacteria in your analysis were one or the other and again do the two tree strategy (point 1).
Point 3 Using large kingdom wide phylogenies from published studies to determine where the most ancestral branch in your taxa.
This is not an exhaustive list about how to use external information to determine the location of the root, but it is very widely recognised.

Directional versus reversible matrices
Directional matrix (e.g. most options in Beast) need to know where the root is, or they will 'guess it': ideally and this information can be feed into the algorithm when it is started. This is important if you are doing molecular dating.
In fact, you can do this even with reversible matrices (e.g. RAxML) but it makes no difference to the answer.

Manual re-rooting
If you have identified the root and the phylogeny you have does not represent that root correctly you can "re-root" the tree using the outgroup identified from point 1 to 3. A program such FigTree are cool for re-rooting. Re-rooting is perfectly cool providing you present the methodology that identifies the root taxa or taxon (singular).

Thanks @M_ for your answer that helps me unclog this issue. I will go through the three approaches to find out the best one in my case. Thanks


text:
How to configure Cromwell to work with slurm?
I have a cluster that is set up with slurm and I seek to use Cromwell to dispatch jobs to the slurm scheduler. I followed the description explained on Cromwells Documentation Website. And created a configuration file called cromwellconfig.conf, however when I run it using the command java -Dconfig.file=/lccs/home/logio9/Configure_Cromwell/cromwellconfig.conf \ -jar cromwell-77.jar run test-workflow.wdl I get an error parsing the config file. I am confused as to where the error is, bellow is the contents of the configuration file:
include required(classpath("application"))

backend {
  default = SLURM
  providers {
    SLURM {
      actor-factory = "cromwell.backend.impl.sfs.config.ConfigBackendLifecycleActorFactory"
      config {
        runtime-attributes = """
        Int runtime_minutes = 600
        Int cpus = 2
        Int requested_memory_mb_per_core = 8000
        String queue = "short"
        String jobname = "testwork"
        String cwd = "/lccs/home/logio9/Configure_Cromwell"
        String err = "errors.output"
        String out = "ouput.out"
        """

        submit = """
            sbatch -J ${jobname} -D ${cwd} -o ${out} -e ${err} -t ${runtime_minutes} -p ${queue} \
            ${"-c " + cpus} \
        --mem-per-cpu ${requested_memory_mb_per_core} \
            --wrap "/bin/bash ${script}"
    """
    kill = "scancel ${job_id}"
        check-alive = "squeue -j ${job_id}"
        job-id-regex = "Submitted batch job (\\d+).*"
      }
    }
  }
}

This is the test wdl script I made, it copies a very large directory into another directory.
version 1.0

workflow Greetings {
    call SayHello {
    }
}

task SayHello {
    command <<<
        cp -R /lccs/home/logio9/Configure_Cromwell/tf1 /lccs/home/logio9/Configure_Cromwell/tf2

    >>>
    output {
        Array[String] out = read_lines(stdout())
    }
}

If anyone has had experience with this, please do share a bit on how exactly this configuration works.

Could you provide what is exactly the error that you are getting?

It says: java.lang.RuntimeException: Error parsing generated wdl: and
Caused by: wdl.draft2.parser.WdlParser$SyntaxError: ERROR: Sibling nodes have conflicting names:
Declaration defined here (line 5, col 1):
String cwd
^

Cromwell already defines the cwd, err and out runtime attributes by itself (it also defines a job_name and script). So when you define these variables in your config it causes a name conflict, since they already exist.
Cromwell will make a directory inside of the "cromwell-executions" directory (created when you run cromwell) for each job it runs. It will write the script, various stdout and stderr files,  and the rc file to this directory. The intend is also for the job to run in that directory. Hence cromwell provides these runtime attributes to point to the appropriate locations within this directory so they can be used in your custom submission command.
You can try just removing these variables in the runtime-attributes section, while keeping them in the submit section.

Thank you for your explanation, your solution works. The only issue that I am having is when I run Cromwell with the command above, Cromwell still seems to be running. If I do a CTRL + C, it stops Cromwell and quits the slurm job. Is there a way to simply dispatch the job to slurm and have Cromwell return/exit? Or perhaps have Cromwell run in the background?

@LogiosJames Cromwell is designed to run multiple consecutive jobs which rely on each other's outputs. So, it needs to keep running in order to check up on the status of the jobs it has dispatched and then submit the next jobs when their inputs have become available. You could just submit the cromwell process to the cluster in the traditional way, though, and have it run on an execution node (eg. sbatch [...] --wrap cromwell run worklow.wdl). Or use something like byobu/tmux and run cromwell in one window and do your other work in another.

Ok, good to note. I normally access and work with Cromwell on my cluster via an ssh connection, if I am running my workflow as I described above and my ssh connection to my cluster closes, would Cromwell continue executing, or would it exit/stop upon the connection closing?

@LogiosJames If you just run cromwell directly in your ssh session it will stop when you disconnect. If you submit cromwell to an execution node, it will keep running on that node. If you open a tmux or byobu session, those sessions will stay alive when you disconnect and you can reconnect with those sessions later using tmux attach-session -t 0 or byobu respectively.


text:
Phylogenetic tree building from proGenomes database for shotgun metagenomics
For some weeks I'm fighting with an issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account phylogenetic distance.
I’m trying to create a phylogenetic tree from the representative genomes of proGenomes v.2.1 database. However, I do not have a representative genome for each taxa id. In contrast, I have a multifasta files for more than 100 contigs for each taxa ID (file freeze12.contigs.representatives.fasta.gz).
For each taxaID (the first numbers in each header of the multifasta file is NCBI taxon ID), I would like to bin the contigs into the biggest scaffold possible, and then perform an alignment (through MUSCLE as examples) followed by RaXML-NG to construct a Newick tree.
However, I cannot find any software for binning contigs with only the contigs as inputs. MetaBAT needs bam files from shotgun read mappings to the contigs in question (not present in the proGenomes database). I found another software (Scaffold_builder), but it needs a representative reference genome. And in microbiology I have more than 100 complete genomes for each NCBI tax id...
For this reason, I ask here an open question if exists some software that can perform of contigs taking into account only a multifasta file of contigs?

I would not recommend the approach that you describe. I am not sure that I completely understand it. For example, metagenomic binning by e.g. MetaBAT is not a scaffolding operation, which returns a single sequence. Rather, they are methods to separate a pool of contigs into groups with similar sequence characteristics indicating they come from the same genome.
It sounds like you already have a group of contigs from each genome, you just want them to all be part of a single sequence in a FASTA file, as would be suggested by your interest in the reference-guided scaffolding approach scaffold_builder.
I would propose that you instead use a tool like GTDBtk to place each of your genomes (each represented by a multifasta of contigs, which is the expected input) into an existing taxonomic tree, and then use their placements in that tree as a model of the taxonomy of your genomes, rather than trying to infer it de novo. I assure you that it is a high quality phylogeny that no one will criticize, and it avoids the possibility of introducing errors into your phylogeny through problematic scaffolding operations.
I would suggest reading GTDBtk documentation for more information.
A narrower answer to the original question
I realized that there is an answer to the question of binning contigs without any additional information. I do not recommend it as I think it would give you positively wrong answers because you are starting with your genomes already binned, and it would not provide you with scaffolds, which I think is what you actually want. But there is a way to perform this task.
You can compute tetranucleotide frequencies from contigs and use them to cluster arbitrary sets of contigs into genomes or other groupings based on similarities. This has a significant rate of error, but it does more or less work. It's an older approach.

Sorry for the delay Maximilian Press. Thanks for your comments. I will try to use the GTDBTk way. I really appreciate your knowledge. I will try to go ahead and overcome this situation. Many thanks.

Thank you another time @Maximilian Press. Yes, GTDB-Tk it's a robust software but although the merging of archaea and bacteria taxa is still a current research topic, I have succeed with another tool that permits me obtaining an unique tree (including archea and bacteria).
If I work with GTDB-Tk I've obtained one tree for archaea and other for bacteria. Then to use it in phyloseq I must work with two independent phyloseq objects which is nice also but my purpose of building a tree from shotgun is only to permits me to calculate diversities based on phylogeny.
As alternative that works it was the following:
to construct a phylogenetic tree from Kaiju I used the 40 marker proteins that is based on the proGenomes database used by kaiju (https://progenomes.embl.de/data/proGenomes2.1_markerGenes.tar.gz) and then I've used the ETE-Toolkit(http://etetoolkit.org/) software to construct the tree in Newick format that I can integrate into a phyloseq object.
I post here if someone can help my strategy. It has been no easy (unless for a dummie as me) to construct a phylogenetic tree from Kaiju's taxonomic classifier.


text:
ATAC peak anatations
I am new in ATAC analysis.
I have a problem about how to transfer ATAC peak to gene symbol.
Here is the ATAC-seq: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM3271041
The GSM3271041_ATAC_sciCAR_A549_peak.txt.gz file contains five columns : (id, peak, chr, start, end). What I am doing is to use bedtools to transfer locations in chromosome to get gene symbols.
Here is an example from GSM3271041_ATAC_sciCAR_A549_peak.txt.gz
id,peak,chr,start,end
9,1-68591-69091,1,68591,69091
10,1-91105-91605,1,91105,91605

I use bedtools to get the gene symbol of id=9 , that is OR4F5.
I want to use these ''symbol-transfered'' ATAC-seq to integrate with RNA-seq which has gene symbols alreadly.
Am I right to do so ?

ATAC-seq data is fundamentally different from RNA-seq data because they measure different things.
RNA-seq data is typically mapped from raw sequencing reads to annotated gene transcripts.
This has a nice 1:1 correspondence because RNA is transcribed from genes.
ATAC-seq data, on the other hand, is used to measure accessible regions of the genome.
But accessible regions don't come from genes the same way that RNA comes from genes.
Accessible regions are scattered throughout the genome and most lie in intergenic regions where there is no annotation for a gene.
Because accessible regions don't have the same biological connection to genes that RNA transcripts do, you need to be more precise in how you want to pair up accessible regions and genes.
Are you looking for cis-regulatory elements (CREs) that regulate nearby genes?
Those will often be found within accessible regions, so you'd need additional data that connects an individual regulatory element to that gene.
Are you looking to create a catalogue of accessible regions within a certain cell type? Then you might want to look at ENCODE's candidate CREs database to see how your ATAC-seq peaks compares.
From there, you'll be able to find additional information about those peaks, such as whether other functional genomics data identifies these regions as active regulatory elements, which gene(s) these elements regulate, and more.
There are different options depending on what exactly you're trying to do, but there is no one way to transfer ATAC-seq peaks to gene IDs because there's no biological question there.

Thanks for your answer!
I think I want to map the peak in scATAC-seq to nearby genes.
I found this.
It's called peak annotation.

Annotating each with with the nearest gene(s) is a more specific question, that's great. Elsewhere on 10X's site you can see a description of their annotation process. Essentially, if you get a list of genes in a BED file, you can use bedtools closest with your peak list to find the nearest genes. The will usually work as a first pass for most applications, but like I said above, it will highly depend on your specific biological question.


text:
"Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification" when running HeatmapAnnotation() (from ComplexHeatmap)
This question was also asked on github
I'm trying to run the function HeatmapAnnotation (a function in the ComplexHeatmap package). When I run it, I get the error "Error in col2rgb(colors, alpha = TRUE) : invalid RGB specification"
col2rgb() is not among the libraries I've loaded and it doesn't seem to be a method contained in HeatmapAnnotation.
Here's the code block for the section:
df <- data.frame(labid= combmetadata$Sample_Name,
             Diagnosis= combmetadata$Diagnosis,
             Tissue = combmetadata$Tissue,
             Tumor.Normal = combmetadata$T.N,
             study = combmetadata$batch)

 labid <- colnames(bValsselect)
 df <- left_join(as.data.frame(labid),df,by="labid")
 df <- df[,c(2,4,5)]
 df$study[df$study == "1" | df$study == "2"] <- "cgid"

 df$Diagnosis <- factor(df$Diagnosis, levels = c("Uterine endometrial","Ovarian endometrioid",
                                            "Ovarian mucinous",
                                            "Colorectal mucinous",
                                            "Pancreas mucinous",
                                            "Stomach mucinous"))

 ha_column = HeatmapAnnotation(df = df,
                          col = list(Diagnosis = c("Uterine endometrial" = "#DDCC77F",
                                          "Ovarian endometrioid" = "#0F7554",
                                          "Ovarian mucinous" = "#44AA99",
                                          "Colorectal mucinous" = "#882255",
                                          "Pancreas mucinous" = "#AA4499",
                                          "Stomach mucinous" = "#D695D0"),
                            Tumor.Normal = c("N" = "grey",
                                             "T" = "black"),
                            study = c("cgid" = "#332288",
                                      "tcga" = "#BD9AF7")))

hegt <- Heatmap(bValsselect,
            column_split=4,
            name = "Beta Value",
            col = colorRamp2(c(0,0.25,0.5),c("#00FFCC","#FFFFFF","#0099FF")),
            show_row_names = FALSE,
            top_annotation = ha_column,
            show_row_dend = F)

The code runs fine up until it reaches the end of the HeatmapAnnotation() call
I've tried to load col2rgb() and override the alpha setting. Can anyone give me a clue about what's happening and any suggestions on how to correct the problem?  I'm running R 4.2. Thanks

Looking at the documentation, it seems that HeatmapAnnotation expects a vector, whereas you are giving it a list. It is not clear on what acceptable color types are. It might even work with a list, but I'd suggest playing around with what you are supplying as your col parameter.
Note that HeatmapAnnotation has to do color mapping of some kind, the code is a little dense, but it does look like it calls out to SingleAnnotation, which in turn runs ColorRamp2, which does use col2rgb.
It also looks like col2rgb is called in ColorRamp2, which you use in the Heatmap call. How sure are you that your script isn't making it that far? It would help to see error output, possibly.
That said, when I take your isolated ColorRamp2 and run it, it goes ok. Possibly pull out the ColorRamp2 call and run it yourself in your environment to make sure nothing weird is happening.
Try moving the assignment of the colors outside the Heatmap call, occasionally supplying function calls as parameters like that causes issues.

You're right about the chain of calls.  I do know that the code didn't run past the HeatmapAnnotation() function because in the process of debugging it, I could generate the error before I run the Heatmap() subblock.  However, I just realized that the problem is that the first color assignment contains 7 alphanumeric characters instead of 6.  That's what was causing the function to fail.  Thanks for your help though.  It certainly helps me understand the function a little better.

The alphanumeric value for the first color assignment contains 7 alphanumeric characters instead of 6. Correcting this typo fixed the problem.


text:
How to manage memory contraints when analyzing a large number of gene count matrices? I keep running out of RAM with my current pipeline
I have several hundred scRNA-seq count matrices, each from a different sample. For my other dataset containg a few dozen samples, I simply merged everything together into one Seurat object, but that won’t work here as far as I can tell. When I try to merge them in the typical way, I get the error:
Error in .cbind2Csp(x, y) :
  Cholmod error 'problem too large' at file ../Core/cholmod_sparse.c, line 89

This basically means I ran out of RAM.  I even ran it on the cluster with 1 TB of RAM, and I got the same error.  This is unsurprising, as memory requirements increase exponentially in relation to matrix size.  It seems that this happens when the merge() function attempts to convert all the sparse matrices to dense matrices.  Here is my code:
seurat_object_list <- "path/to/matrices" |> list.dirs(recursive = FALSE) |> lapply(\(matrix_folder) {
    matrix_folder |> Seurat::Read10X() |> Seurat::CreateSeuratObject()
})
seurat <- merge(seurat_object_list[[1]], seurat_object_list[2 : length(seurat_object_list)])

The Seurat object list is created with no issue, but it always fails on the merge.  I tested this code with a smaller subset of the matrices, and it works just fine.  It seems that I just have too many cells to process in this manner.
At this point, It looks like I have a few options.

I could hack together some way to force the matrices into one Seurat object.  But I don't know if that would even be useful, since I would likely keep bumping my head on the same memory contraints whenever I try to do anything with the combined object, such as running SCTransform() or FindMarkers().
I could divide the samples into groups that are small enough to merge, and process them in the normal way.  But this imposes limitations on the types of comparisons/analyses I could perform.
I could keep the matrices separate and analyze them iteratively.  This would solve the memory problems, but I haven't found any info online about anyone analyzing count matrices this way.  I imagine that I would lose a lot of out-of-the-box functionality that comes with the standard single-cell R libraries.  For example, I'm not sure how I would normalize all the gene counts against each other as separate objects.  This option seems like it might require a lot of re-inventing the wheel and a lot of code.

I like option 3 the best, but I'm not sure where to begin.  I am still pretty new to bioninformatics, so I would love to hear some input from those who are more experienced.  How do you deal with memory constraints when dealing with a large number of cells/samples?  Are there any standard ways to deal with this problem?  Is there some tutorial/doc online that I failed to find?  Any advice you could offer is greatly appreciated.

I get the same error but not running out of memory. I don't think this is a memory problem. It's some internal limitation of dgCMatrix object.

Scarf has been developed to specifically solve this issue. Scarf handles  large-scale scRNA-Seq/scATAC-Seq/CITE-Seq datasets that will otherwise not fit into memory. Unfortunately, this is not available in R yet.
Here is the Github link:
https://github.com/parashardhapola/scarf
Documentation:
https://scarf.readthedocs.io/
Paper:
https://www.nature.com/articles/s41467-022-32097-3
Please feel free to create an issue on Github for you use case (if you cant already find a solution in the documentation).
/PD


text:
Creating msk file for Illumina Array Analysis CLI Copy number train
recently I'm researching on how to detect CNV from Illumina microarray data. I found that Illumina has Array Analysis CLI software for this task, and I need to train the CN model myself, since the official .dat file is no where to be found.
However, this requires another file: the manifest file (.bpm) must be accompanied by a "mask file" (.msk), which I cannot find online.
Could you please help me on this issue please? And could you suggest some CNV detector for Illumina microarray? Thank you.


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

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

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


text:
Does the number of RNA reads per cell obtained from the 10X scRNA experiment depend on amount of mRNA in given cell?
As we know, the amount of RNA reads per cell obtained from 10X scRNA experiment vary between cells. I wonder if this is effect of technical issues or does the number of RNA reads per cell obtained from the 10X scRNA experiment depend on the amount of mRNA in that cell? For the latter to be the case, we would have to read almost all the mRNA fragments in a cell, which I don't think is likely.
Thank you for your replies in advance!

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

I am inclined to agree with you that the reads-per-cell is dependent on the cell type and cell cycle state.  In the datasets I have analyzed, I observed that immune cells tend to have far fewer reads than fibroblasts, which have fewer reads than larger epithelial cell types.
EDIT: I have also observed that macrophage populations generally have more RNA than lymphocytes. Granulocytes have even less if you even capture them in your sample.

The number of reads per cell for any given experiment is largely technical. It depends on sequencing depth. 10X recommend 20K reads per cell for the v3 chemistry (and 50K reads per cell for v2)

Thank you for your answer. I'm aware of that. What I wanted to ask is:
If in my data from one 10X experiment and from 1 sequencing I have macrophages cluster that have 100k reads per cell and lymphocytes cluster with 33k reads per cell, can we assume that the amount of mRNA in these lymphocytes cluster was 3times smaller (in living cells, befor 10X exp) than in given macrophages cluster?


text:
How to build advanced boxplots from datasets which columns contain only numbers?
I have got a fold-change dataset and I simply would like to plot something like this:

The problem is that I have only numbers in my columns and when I use a boxplot() function I can plot only boxes but not dots, jitter, etc.
I use this code to build a df and then proceed to a boxplot
A <- c(-1, 25, 3, 4, 50)
B <- c(23, -41, 32, 68, 26)
C <- c(63, 71, 22, -8, 2)
D <- c(16, -49, 2, 56, 6)
df <- data.frame(A, B, C, D)
print (df)
boxplot(df)

How can I plot something that looks like the above picture? Thanks

You want ggplot2, specifically the geom_boxplot and geom_dotplot functions. Check out the cheat sheet (warning: PDF download) to see the full functionality of the package.

Thank you. The problem is that I don't have geom data - if you check the dataset I have proposed - there is no "y" data, I have only 4 columns (A,B,C,D) and rows with numbers. It is ok for default boxplot() while I cannot add geom, can you specify what should the ggplot2 script look like for the above-mentioned dataset? Thanks

To coerce your data into a form interpretable by ggplot2, you can use the gather function from tidyr.
gather has been superseded by pivot_longer.
library(tidyverse)
df <- tibble(A = c(-1, 25, 3, 4, 50),
             B = c(23, -41, 32, 68, 26),
             C = c(63, 71, 22, -8, 2),
             D = c(16, -49, 2, 56, 6))

tidyr::pivot_longer(data = df, cols = c(A,B,C,D)) |>
  ggplot(mapping = aes(x = name, y = value)) +
  geom_boxplot() +
  geom_dotplot(binaxis = "y", stackdir = "center") +
  theme_classic()


text:
Cytoscape for coloring the nodes in the network
I need to color different categories of the proteins in the interaction network automatically. Is there any way to color each category with a different color in cytoscape? or R?
My input is 2 column interaction protein ID. I have different categories of proteins in CSV format separately only with protein ID.

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

In protein -protein interaction network, I need to color the proteins of the different categories with different colors. I have four categories of proteins.

https://stackoverflow.com/questions/15999877/correctly-color-vertices-in-r-igraph, https://manual.cytoscape.org/en/stable/Styles.html


text:
What does this accession NCBI code mean: 6MWN_B?
According to this article, accession codes should consist from a combination of uppercase letters following a combination of digits. If this is a RefSeq, it can have a prefix as a combination of uppercase letters with underscore.
But this accession code 6MWN_B does not consist with all aforementioned. What is it? I could not find an answer via Google.

In general, an accession is a value that points to (provides access to) a single database entry. Different databases (GenBank, RefSeq, Ensembl, SwissProt, and so on) have different conventions for accessions.
The accession in question points to an RNA sequence associated with a Protein Data Bank (PDB) record: "Crystal structure of hepatitis A virus IRES domain V in complex with Fab HAVx". In the PDB record, you can click the "sequence" tab and select chain "B" from the drop-down, and you will see an RNA sequence (ACGU) matching the DNA sequence (ACGT) in the NCBI record to which you linked.
I'm not sure what the rules or conventions for PDB accessions are, but they are clearly different from RefSeq accessions.

Thank you! I am just surprised that this lies in nuccore & did not get a new local accession code.


text:
Where to download a table with ICD-9-CM codes?
I am looking for a simple table with all ICD9 codes (For example: Diagnosis Code 414.0 -> Coronary atherosclerosis). I can only find weird PDF or other kinds of text documents. Or web forms, where I have to post a certain code to get the description.
Here is what I have found so far.
http://www.icd9data.com/2015/Volume1/default.html
https://www.aapc.com/codes/icd9-codes-vol3-range
https://www.cms.gov/Medicare/Coding/ICD9ProviderDiagnosticCodes/codes
https://www2.gov.bc.ca/gov/content/health/practitioner-professional-resources/msp/physicians/diagnostic-code-descriptions-icd-9
https://ftp.cdc.gov/pub/Health_Statistics/NCHS/Publications/

I wish to access icd-9cm-2015 as an excel file. Any information or help would be welcome.

Hi @monapanahifar and welcome to the site. You might think about opening a separate question with a link to this one. Please note I did edit your comment a little. Please check that is what was intended. If not please drop a note in chat.

The zip files found within https://www.cms.gov/Medicare/Coding/ICD9ProviderDiagnosticCodes/codes contain text files and Excel spreadsheets that you can open with Excel, R or Python depending on how you want to use this data.

Thanks, I saw those. However, they seem to have a different format. I will add an example.


text:
Extract SNPs from multiple sequence alignment
I’m wondering if folks have recommendations for tools/scripts to extract SNP sites from a multiple sequence alignment of consensus bacterial genomes? Specifically, I am interested in a multiple sequence alignment containing variant (SNP) positions only, to be used as input for phylogeny-building programs including RAxML and BEAST.
I was previously using https://github.com/sanger-pathogens/snp-sites, but this is not being maintained and I’m running into errors that it seems many others are having (with regards to sequences being different lengths, although I have confirmed they are the same length).
Any recommendations would be great, thank you!

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

Are you looking for identity, or are you variant calling, or finally perhaps looking a consensus sequence?

Yes, apologies for the lack of clarity. I have aligned consensus sequences, I'm hoping to extract SNP sites only into an alignment of SNP positions. Thank you!


text:
What is the difference between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus for converting to gene symbols?
I am interested in the differences between AnnotationDbi/org.Mm.eg.db and biomaRt/Mus.musculus. They yield the same results as seen in the code below.
library(edgeR)
library(biomaRt)
library(Mus.musculus)
library(AnnotationDbi)
library(org.Mm.eg.db)

url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb")
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
           "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt",
           "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt",
           "GSM1545545_JMS9-P8c.txt")

x <- readDGE(files, columns=c(1,3))
geneid <- rownames(x$counts)

genes <- AnnotationDbi::select(org.Mm.eg.db, keys=geneid, columns="SYMBOL", keytype="ENTREZID")
dim(genes)
genes2 <- biomaRt::select(Mus.musculus, keys = geneid, columns = c("SYMBOL"), keytype = "ENTREZID")
dim(genes2)

Do they use different databases? Is one preferred or recommended over the other?

AnnotationDBi uses AnnotationDb objects which is the virtual base class for all annotation packages. It contain a database connection and is meant to be the parent for a set of classes in the Bioconductor annotation packages. Examples of classes derived from the AnnotationDb object include (but are not limited to): ChipDb, OrgDb, GODb, OrthologyDb and ReactomeDb. The database from one of these classes is used to extract annotations. org.Mm.eg.db belongs to OrgDb class.
Mus.musculus  represents a collection of annotation packages that can be used as a single object. It combines different information about same organisms (Mus musculus) into one place. It also uses org.Mm.eg.db as a dependency. So, both uses same database to extract Gene symbols. That is the reason they are producing same results. There is helpful documentation available on Bioconductor.


text:
chip seq biological replicate merge after peak-calling
I have 3 biological replicates for KO condition. What I have done so far

Peak calling
The removed blacklisted region

Now Im bit confused about the peak merging stage. Here Im trying to use the bedops tool.
So what is the best option in my case when I want to make a combined set of peaks that is present in all the three replicates what option do i have to use?
Which option to use from bedops
 Everything (-u, –everything) or `Merge (-m, –merge)`

consider calling peaks with Genrich which will handle you replicates quite nicely.

although I haven't used it but will give it try


text:
Violin plots appear as vertical lines
I want to perform a violin plot from a Seurat object and I rather have vertical lines. How can I stop this and rather have a standard Violin plot?
code
Seurat::VlnPlot(object, features=c("LTA", "LTB"), group.by = "group",
    pt.size = 0)

Why have you set the point size to zero? I expect that would stop the points from displaying on the plot.

I believe the reason is simply LTA and LTB not being expressed (or their expression not being detected due to technical limitations). If you would set the pt.size argument to something else than 0, you would see a big blob of dots at point 0.0 (depending on the number of cells/data points).


text:
Association analysis for methylation beta value and phenotypes in R
I am trying to learn how to do association analysis of methylation array data with numerical trait like body weight, eGFR.
As I have seen some of the publications, methylation data association with eGFR has been showed through regression analysis but I am not sure how to run regression for so many CpG probes in R.
Can you have look on this below publication and could you help with R code for regression analysis?
https://www.nature.com/articles/s41467-017-01297-7.pdf
I do have methylation beta value from RnBeads but I am not sure how to read methylation data (rows/column-CpG/SampleID) and then how to read phenotype data file and finally how to make single sheet to run regression analysis?

Have a look at this method https://academic.oup.com/bioinformatics/article/34/15/2657/4939328

Hi, I found a way to do methylation data association analysis with eGFR as shown here (https://support.bioconductor.org/p/93868/), but not sure how make input files; like in methylation data file; do I take CpG probe as rows or column? and I am also not sure how to mention covarites in modal.matrix for eGFR association? Many thanks

Answer from @user4939328, converted from comment. Please improve this answer if you can:
Have a look at this method:
https://academic.oup.com/bioinformatics/article/34/15/2657/4939328
[Wondering why this was converted from a comment? See the StackExchange rules about comments here]


text:
What is the best way to perform batch ajustment between two experiments?
I have duplicate drug efficacy experiments performed six months apart. I would like to know what is the right way to perform statistical analysis to understand fold changes after treatment with a drug.
Say I have
ctrl1       tube1       ctrl2       tube2
7,732637722 10,71639307 5,454318175 8,052899428
6,69630905  8,762916793 7,441582994 9,855106201
10,35624396 11,73716043 8,976524001 11,0975511
4,29150392  4,848632578 5,901132154 7,918630661

and the same experiment but 6 months after (numbers are slightly different)
Could I take the mean between samples divided by two can work out? If not - why?
I would like to have statistically significant values for futher analysis.
Thanks.

What do you want to show? Taking the mean is rarely a good move except in the final analysis

I want to adjust values between the same experiments, to be able to analyze them via t-test for example. It's just these two batches are 6-months different.

The answer is simply you use both sets of values in a one-way paired ANOVA. Note that certain criteria need to be in place for ANOVA, notably the Fmax test.
Looking at the data 2-way ANOVA might be needed and that is more complicated, or else 2 one-way ANOVAs: I don't understand 'ctl2' and 'tube2' data. Is this a single or multiple drugs?
Wilcoxon signed-rank test is the first choice test because there are no prior assumptions required, which replaces both the T-test and the 1-way ANOVA.
Whatever you do, taking the mean and performing a T-test is not a good idea. You asked why ... because the variation between the values is critical to the significance of the test and is simply demonstrated:
Tube 1 value
tube1 time X = 11
tube1 time y = 9
tube 1 mean = 10

Tube 2 value
tube 2 time x = 19
tube 2 time y = 1
tube 2 mean = 10

Thus both would appear to have no change - they are both 10 - and yet the variation between them means they observe very different behaviours. Thats why.


text:
how to create "sample file" for the qAlign() function after trimming the reads in R
I'm an absolute beginner trying to solve this question "Align the trimmed and untrimmed reads using QuasR and plot alignment statistics, did the trimming improve alignments?"
I did trim the reads using the preprocessReads() function, now I want to align the output but the qAlign() only takes a samples file and genome file. however, my trimmed reads are two separate files how can I align them.
I asked google but I didn't find the answer

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

Please consider the answer below and whether you wish to upvote/accept.

The QuasR vignette tells you how to run it. The samples file should contain a list of your samples and fastq files.


text:
How to use EPAM Cloud Pipelines?
I am searching for an easy-to-use GUI management tool to handle my WDL workflows, I found EPAM Cloud Pipelines. Unfortunately, it uses cloud services such as AWS, GCP, etc. but I wanted to use my own backend instead, thus, I attempted to see if it would be possible to configure it to use my backend, however, there seems to be little to no instructions on how to install it. If anyone has any experience using this, can you please provide a description? Or, if you use or know of another GUI tool similar to EPAM Cloud Pipelines that can help manage and build WDL workflows, please share some information on that aswell.


text:
How can you put a bash array into a wdl variable?
Given two separate wdl arrays, my command section will do some modifications to each array and place the output into two separate bash arrays. Finally, I want to output the two separate bash arrays in the output section of the wdl task to be used in other tasks within the workflow.  The below code is an example of what I would like to do:
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        bashX = ("Truea" Falsea")
        bashY = ("Gooda" "Fooda")
    >>>
    output {
        Array[String] finX = bashX
        Array[String] finY = bashY
    }
}

The above example produces an error. It seems that I might have to convert the bash arrays back into wdl arrays in order to declare them as outputs in the output section, however, I am not sure how to do this. May someone help me with this? Moreover, a general explanation of how you can capture the output of specific things from the command section into the output section would be helpful aswell.

When a WDL task is run, the execution engine takes the command block and evaluates it as a template, rendering any wdl expressions it contains into the final script that is executed, BEFORE it is actually executed.
The actual execution environment is generally isolated from the execution engine or from other tasks. Information is passed between tasks by reading from or passing around file based outputs and not from accessing the bash environment directly, which the engine does not have access to.
So, in order to achieve what you are trying to do, you first need to get the values into a file format, and then you need to read those values using an engine function provided by WDL. Generally, all outputs need to be in a file if you want the engine to be able to interact with them. Please note that stdout and stderr are considered "files" and are accessible with stdout() and stderr() engine functions.
Option 1: Construct json files and use read_json
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:
        echo '["Truea","Falsea"]' > x.json
        echo '["Gooda","Fooda"]' > y.json
    >>>
    output {
        Array[String] finX = read_json("x.json")
        Array[String] finY = read_json("y.json")
    }
}

Option 2: use read_lines
task AddA {
    input {
        Array[String] wdllistX = [ "True", "False"]
        Array[String] wdllistY = [ "Good", "Food"]
    }

    command <<<
        ###We add the letter 'a' at the end of each string in each wdl list and finally get the
      bash arrays:

        echo "Truea" >> x.txt
        echo "Falsea" >> x.txt
        echo "Gooda" >> y.txt
        echo "Fooda" >> y.txt
    >>>
    output {
        Array[String] finX = read_lines("x.txt")
        Array[String] finY = read_lines("y.txt")
    }
}


text:
Drug indication data for only cancer phenotype from Chembl
I am interested in getting the drug indication data for all drugs in the database using sqlite.
I used following step to get this
import pandas as pd
import sqlite3
con = sqlite3.connect("chembl_30.db")
df = pd.read_sql_query("SELECT * FROM drug_indication", con)

With this I am getting a table with columns efo_id and efo_term for all disease.
drugind_id  record_id   molregno    max_phase_for_ind   mesh_id mesh_heading    efo_id  efo_term
0   22606   1390869 675774  2   D045743 Scleroderma, Diffuse    EFO:0000404 diffuse scleroderma
1   22607   1390869 675774  4   D001172 Arthritis, Rheumatoid   EFO:0000685 rheumatoid arthritis
2   22608   1390869 675774  2   D008180 Lupus Erythematosus, Systemic   EFO:0002690 systemic lupus erythematosus
3   22609   1343491 675576  3   D009203 Myocardial Infarction   EFO:0000612 myocardial infarction
4   22610   1344515 438435  4   D011471 Prostatic Neoplasms EFO:0001663 prostate carcinoma

Is there any way I can I filter EFO's based on cancer, so that I get drug indication data (max_phase) for all chembl ids used in cancer?
some important links
https://www.ebi.ac.uk/chembl/db_schema
https://www.ebi.ac.uk/chembl/compound_report_card/CHEMBL941/
Thank you


text:
VCF or BAM file for raw data of gene test?
My friend has a VUS (Polymicrogyria with or without VEDS) mutation that was found in her whole-exome sequencing with respect to the phenotype given that time. At present, her doctors have more information on the complexity and are looking for more explainable mutation present in herself, if any. We asked the company and they're willing to give the raw data.
Which file format should we opt for? They've given the option of BAM & VCF.

First of all - I would advise you against concluding anything out of your own analysis of the data, you need to consult someone with appropriate training!
Can you get both? The .bam file is a format that stores information about sequencing reads (usually mapped, but not always). You can even extract the original sequencing information from bam files. This is practically the rawest information you can have. The .vcf file is a variant file that gives you some information about inferred variants, but there are plenty of bits of information that will be missing. The basic information regarding sequencing depth and quality of variant calling is still included and therefore most likely sufficient for your physician to see any other variants potentially relevant to your friend's condition.
I will repeat once again - don't conclude anything on your own. It takes years to gain expertise in analysing genemics data and there is very high risk of misinterpretations if not analysed with care.


text:
Python; Question about the append (modifying lists): why modifying a list with append changes others?
I have a list named dict_eg1 whose values are lists.
dict_eg1 ={
    "A":[1,2,3],
    "B":[11,22,33],}
val1=dict_eg1["A"]
val1.append (100)#modify the list by using append.
print (dict_eg1)

{'A': [1, 2, 3, 100], 'B': [11, 22, 33]}
In this case, dict_eg1 changed.
However,
dict_eg1 ={
    "A":[1,2,3],
    "B":[11,22,33],}
val1=dict_eg1["A"]
val1 = [1, 2, 3, 100]#modify the list without using append.
dict_eg1

{'A': [1, 2, 3], 'B': [11, 22, 33]}
In this case, dict_eg1 did not change.
Can anyone tell me why these two results were different??

I’m voting to close this question because it is a pure python question and not related to bioinformatics.

To sum up, .append() modifies the list in-place, while assignment with = creates a new list entirely

Welcome to Bioinformatics.
First of all, it should not be asked in this community since it is purely a question about python programming. I suggest that next time, please head to StackOverflow as you will likely receive more help from there.
You also look new to programming. This is more of a logical question. Let's break it up to understand what is happening in your first example:
# you define your dictionary containing two keys along with their values.
dict_eg1 = {
             "A":[1,2,3],    # Key 'A' with a value as a list [1,2,3]
             "B":[11,22,33]  # Key 'B' with a value as a list [11,22,33]
}

# you now define a new variable 'val1' which is
# the value of key 'A' and that is the list [1,2,3]
val1=dict_eg1["A"]

# as you have mentioned, doing the below command
# will append the value to the list of key 'A'
# in other words, you can also write it like so
# dict_eg1["A"].append(100)
# that's why when you make a change to 'val1'
# the value of the key 'A' of 'dict_eg1' will change
val1.append (100)

As for the second example (I will skip the same steps):
# ...

# you store the value of key 'A' in 'val1'
val1=dict_eg1["A"]

# however, in the below command, you actually
# !redefine! the value of 'val1' to a new list
# that is totally independent of your 'dict_eg1'
# In this case, 'dict_eg1' is not affected by
# any changes made to 'val1'
val1 = [1, 2, 3, 100]

Just to sum it up, when you use val1=dict_eg1["A"], imagine that you link the value of your key A to val1 and any changes made to val1 will affect dict_eg1["A"]. However, redefining val1 will disconnect that link.

Thank you. I am sorry that I posted this question in the inappropriate community.


text:
Different results of DESeq2 in different ways of comparison
I have different treatments (infected and control ones). When using DESeq2, I have different results in two ways of comparison.
first way:
factor: treatment, factor level1: Moribund Infected, factor level2: 5week Infected, factor level3: Moribund Control, factor level4: 5week Control.
I have 6 result files (comparisons among factor levels)
second way: (I compared just 2 groups to each other)
factor: treatment, factor level1: Moribund Infected, factor level2: 5week Infected.
factor: treatment, factor level1: Moribund Infected, factor level2: Moribund Control.
factor: treatment, factor level1: Moribund Infected, factor level2:the  5week Control.
factor: treatment, factor level1: 5week Infected, factor level2: Moribund Control.
factor: treatment, factor level1: 5week Infected, factor level2: 5week Control.
factor: treatment, factor level 1: Moribund Control, factor level2: 5week Control.

Finally, the results were different in these two ways of comparison. I have more significantly expressed genes in most comparisons of first way and even I have few comparisons in the second way that have more significantly expressed genes.
Why and which way is trustable?

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

Could you please show the specific design formulas you used to describe the statistical models in the DESeq2 commands. Showing the full factor table would be helpful as well (e.g. colData in the DESeq2 commands).

Cross-posted:https://support.bioconductor.org/p/9144330/


text:
ggtree：How to keep leading zero ‘0’ in the tiplabel？
I am using ggtree to draw a phylogenetic tree.
Some taxa IDs have a 0 at the beginning.
However, when I plot the tiplabel the leading zero is automatically dropped. How do I retain the leading zero?
Thanks.
The code and tree file and metadata file are attached. I would appreciate any help.
setwd('G:\\Experiments\\Paper\\Torula\\06_iqtree')
dir()

# Load required packages
library(xlsx)
library(ggtree)
library(treeio)
library(ggplot2)

# Read tree
tre <- read.iqtree("concatenate_iqtree2.treefile")
tre_tbl <- fortify(tre)
# Read associated information table
metadf <- read.xlsx('G:\\Experiments\\Paper\\Torula\\00_data\\Torula_taxa.xlsx', sheetIndex = 1)

# Plot tree
ggtree(tre) %<+% metadf +
    # our own taxa darkblue
    # type strain bold
    geom_tiplab(aes(subset = (New == 1 & Type == 1),
                    label = paste0('bolditalic(', Genus, ')~bolditalic(', Epithet, ')~bold(', Collection,')~bold(', Number,')')),
                color = 'darkblue', parse = T) +
    geom_tiplab(aes(subset = (New == 0 & Type == 1),
                    label = paste0('bolditalic(', Genus, ')~bolditalic(', Epithet, ')~bold(', Collection,')~bold(', Number,')')),
                color = 'black', parse = T) +
    geom_tiplab(aes(subset = (New == 0 & Type == 0),
                    label = paste0('italic(', Genus, ')~italic(', Epithet, ')', '~', Collection,'~', Number)),
                color = 'black', parse = T) +

    geom_text(aes(label= UFboot), vjust = 0, hjust = 1, size = 3) +
    xlim(NA,0.5)

Some taxa in the figure

I like both the methods and the study. ggplot2 is R's standout package/library and the library ggtree looks great.
Anyway, I would know how to solve this in Python or Perl and I'll use the same logic for R. My suspicion is R is automatically type casting your metadata and therefore has no option but to remove the zero.
Simply convert your integer range into a string, viz.
 Number2 = as.character(Number)

or simply in situ
as.character(Number)

This doesn't make any difference to your output, who cares whether its and integer or string and forces R into retaining the zero, '0'.

Note below that tried and didn't work. A string should be immutable.
The suspicion is the the '0' is being removed at an earlier point in the code, i.e. as soon as it type casts it - zero disappears.
What I would do some print outs of the variable between its import and writing to disk. What I suspect is the conversion to string needs to occur earlier.

Second idea
The other possibility is to wrap the metadata in double/single quotes in the Excel file (or flat file). This way R cannot force it into an integer when it is imported and must retain it as a string. I'd be amazed if R could wriggle out of that.

I already tried that, but it did not work. Thanks anyway.

I've addressed this above @YPCHEN. I'm pretty certain this is a type cast bug and no more.

I quote the metadata in the excel, it works. Thank you very much.

Glad it works @YPCHEN please do remember to 'accept' the answer.


text:
Why is the p-value significance threshold for HLA association tests $5*10^{-8}$?
Typically, a p-value of $p<=5*10^{-8}$ is used for genome-wide association testing, as there are roughly $1/p$ independent LD blocks in the human genome, so it correct for multiple testing across the genome.
HLA assocation tests, such as in this study, only look at a fraction of the genome, something like chr6:28–34 Mb, which is approximately 0.2% of the genome used to conduct a true 'genome-wide' association study. However, they still use the same genome-wide p-value threshold of $p<=5*10^{-8}$.
Is this not too conservative, since there is a much reduced multiple testing burden? Should the threshold not be something like $p/0.002$?

Its a very good question and it is also a good study. What you are saying must be correct.
I thought it was possible they used the correct critical probability in the original submission and a reviewer requested it was changed. However I checked their medRxiv upload here and the authors clearly use this threshold prior their formal submission.

We performed stepwise conditional analysis 479 to identify additional
independent signals by adjusting for the most significant amino acid
480 position in each step until none met the significance threshold (
P = 5 10^-8)

Thus I don't know.
If this paper is central to your area of interest I would use the corrected critical probability threshold and re-assess and extend their conclusions because 5 decimal places is a huge a difference.

As a direct answer to your question, I would argue that the authors, in choosing this threshold, have demonstrated a lack of understanding of their process of analysis. But that answer is not a full answer, and it's a problematic answer because [almost] everyone who does population genetics does similar things. This p-value is not appropriate, not merely because it is at a genome-wide significance level, but because no threshold is appropriate.
To explain my perspective more, thresholding to identify "statistically-significant" results, in general, is not a good approach. Regardless of what threshold is used, p-values only indicate the confidence in a statistic's value (as it relates to fitting a particular statistical model), not importance or the size of the value.
The authors of this paper, as with many other authors who have written similar papers, make an obvious statistical error in the legend of figure 3, where they state: "One-field classical allele HLA-B*57 (P = 9.84 × 10−138) (a) and amino acid position 97 in HLA-B (Pomnibus = 1.86 × 10−184) (b) showed the strongest association signal." This is an error, because they have only reported here the p-values, and p-values have no relevance to the strength or magnitude of a statistic.
Given this, it doesn't make sense to rank by p-value, and it especially doesn't make sense to display p-values on a graph and highlight p-value peaks as indicating association peaks. Those regions are confidence peaks.
In an attempt to provide some constructive light around this, I will answer the obvious follow-up questions are "How should p-values be used?" and "If p-values shouldn't be displayed or ranked, then what should?"
My preference is to use p-values as a loose filter, to exclude obviously incorrect results. In the case of a GWAS, I might, for example, choose to exclude any statistic that has a p-value above 0.1, leaving the remainder as potentially correct. It could possibly be argued that this is still thresholding. It is, in a sense, but importantly I am not saying that all values that pass the threshold are good, just that the ones that are excluded are all bad; I allow for some rubbish to end up in the good pile. I liken it to using a computer to find the things that look like diamonds, so that I don't waste my time inspecting dirt.
After that filtering is done, the p-values should be ignored, and the descriptive statistic used to rank, determine significance, association, or whatever the goal of the analysis is. For GWAS, a commonly-used value is beta, which can be positive or negative depending on the direction of association, and relates to the variant difference between groups.
If a list needs to be cut down, then this descriptive statistic should be used for thresholding, depending on the analysis. For example, it might be appropriate to display the 100 most extreme beta values, or threshold to exclude any beta values that have a magnitude less than 0.01, or display a manhattan-like plot with beta on the $y$ axis, and chromosome position on the $x$ axis.
There shouldn't be set rules, because every situation is different - there's a human element to this. What I'm saying is that it should be the job of the researcher, not the computer, to identify which data are significant. Make the data set small enough that it can be processed by a person, then let them do the rest. I think Hannah Fry explains this quite well in this talk:

Now this is the stuff the algorithms are amazing at: looking for tiny, tiny clues in seemingly completely disconnected data sets.... Just incredibly, incredibly sensitive these algorithms, right? Just absolutely amazing what they can do. And yet remember that sensitivity does not make a perfect algorithm. That's only one half of the equation.... [Humans] are terrible at sensitivity, right? We are really, really bad. We miss things all the time. But, specificity, being specific, that's like our superpower.... And rather than choosing between human and machine, which is kind of the rhetoric that we get so much, why don't we exploit each other's strengths and just create much more of a partnership?.... Because the algorithm never gets tired, so let it trawl through all of that data and just highlight a few key areas of concern. And then the human never misdiagnoses. So they can come in, and just sweep up and have the final say.


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

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

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

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

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

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

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


text:
How to extract reads that map exclusively to a single site with 1 or zero mismatches from BAM files
I generated BAM files (sorted by coordinates) by aligning human RNA reads against the human reference genome using BWA MEM and TopHat2. I would like to know how I can retain only reads that mapped exclusively to a single site in the reference genome with one or zero mismatches for downstream analysis.
For retaining reads mapped to only one site I tried the following two commands, using grep
samtools view -h file1.sorted.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' |
samtools view -b > file1.filtered_unqiue.sorted.bam

or filtering by mapping quality
samtools view -q 10 -b file1.sorted.bam > file1.filtered_unique.sorted.bam

Then I filtered file1.filtered_unique.sorted.bam of either two commands above for mismatches with
samtools view -e '[NM]<=0' -O BAM -o
file1.filtered_unique_mismatches.sorted.bam
file1.filtered_unique.sorted.bam

I am quite unsure about the correctness of the results.
I also came across a one-liner combining both efforts, where the BAM is first intersected with the coordinates of interest and then filtered for mismatches:
samtools view -h file1.sorted.bam chr:start-end | egrep '(^@|(NM:i:(0|1)))' > file1.filtered_unique_mismatches.sorted.bam

However, here I think this would also include multiple mapped reads since I only look for reads mapped in a specified region?
Does anyone have an advice what the best solution is?

Just a small tip: to help you see if your results are correct, you can visualize the resulting bam on a genome browser like IGV.

If you can program in Python, for this kind of task, I would recommend using pysam.


text:
set -euo pipefail in Snakemake wrapper
I have snakemake script doing fastqc with a wrapper but it errors out saying
Command 'set -euo pipefail;  /home/mi/ppp/miniconda3/envs/meta/bin/python3.9 /workflow/.snakemake/scripts/tmpsl695gpr.wrapper.py' returned non-zero exit status 1.
  File "/workflow/rules/fastqc.smk", line 32, in __rule_fastq1
  File "/home/mi/ppp/miniconda3/envs/meta/lib/python3.9/concurrent/futures/thread.py", line 58, in run

my rule looks like this
rule fastq1:
    input:
        "../resources/samples/{sample}_1.fq.gz"
    output:
        html = "../results/html/{sample}_1.fq.html",
        zip = "../results/zip/{sample}_1.fq.zip"
    conda:
        "../envs/meta.yaml"
    log:
        "../results/logs/fastq_raw/fastq_raw_{sample}.log"
    threads:
        2
    wrapper:
        "v1.5.0/bio/fastqc"

I have a feeling this used to work before and my partner says it works fine on his machine but I can't find the issue.
I'm also doing this on a clean conda/mamba install and env.
EDIT:
Snakemake version is 7.7.0
I "solved" this by installing fastqc with conda... I mean am I stupid? Is it a requirement for a wrapper to also have the package installed? I don't think this was mentioned on the Snakemake website
(Im not sure wht mr communuty bot is complaining about the problem is the error message and me not being able to run the code..)
Thanks for reading.

Please include Snakemake version in question body.

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

Note that set -euo pipefail is something Snakemake adds by default to your shell commands. The fact that it appears in the error message is usually not relevant. The error likely has to do with something else in the command that the wrapper calls.


text:
Different colors in UMAP
I have integrated top ten TCR clones into Seurat object something like below

You are seeing some clones are expressed in more than one cells
So in Seurat I have something like
> head(pbmc@meta.data$highlight)
[1] NA                              "CAGVDSNYQLIW_CASSPVAGGSYNEQFF" NA
[4] NA                              NA                              NA

Based on the number of cells in which each clone is expressed, I want to rank clones and says how percentage of cells expressed each clone these figures below

Can you help me with this please?