text:
Confusing result from Sanger Imputation Service (Eagle v2.4 for phasing, PBWT v3.1 for imputation)
I ran my 23andMe data through the Sanger Imputation Service that uses Eagle v2.4 for phasing, PBWT v3.1 for imputation.
However, some aspects of the results are very confusing to me. Perhaps see #3 (the most confusing) first.
Point 1

There are 62,623 rsIDs for chromosome 3 in my 23andMe file and 2,821,894 SNPs (actually, unique rsIDs) after imputation. However, only 60,144 of the original 62,623 rsIDs are in that ~2.8 million. Where did the other 2,479 go? Here are some examples of missing rsIDs:
rs1516332,
rs12374025,
rs1170695,
rs2731343 and
rs63749864.

EDIT Point 1
Looking at rs1516332, the position is 3:105793, which does have a result in the file, but it doesn't have an rsID:
3    105793    .    C    T    .    PASS    RefPanelAF=0.17878;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

Looking at rs63749864 (3:13564048) it's missing by rsID and position.
EDIT Point 1 (again)
As described here: https://github.com/richarddurbin/pbwt/issues/51

"... of the 2,479 'missing' SNPs, 1,806 can be 'found' [by position] leaving 673 'anchors' [or rsIDs called by 23andMe] missing from chromosome 3."

Point 2

Some imputed rsIDs appear more than once in the results. How should I interpret that? e.g.

3   64071   rs116059446 G   T   .   PASS    RefPanelAF=0.000369572;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   64071   rs116059446 G   A   .   PASS    RefPanelAF=9.2393e-05;AN=2;AC=0;INFO=1  GT:ADS:DS:GP    0|0:0,0:0:1,0,0

EDIT Point 2
From the answer below, this is a variant with more than two alleles: rs116059446
Seems the format of this VCF file only allows one alt per line, hence the two alts are split over two lines.
See my question here: Is it valid VCF not to 'squash' positions with more than one ALT allele?
And the answer from Durbo here: https://github.com/richarddurbin/pbwt/issues/52#issuecomment-1062951105
In the above case, the genotype is easy to interpret (REF/REF in both lines).
However, I'm now finding examples that I just can't interpret:
3   255395  rs331869    C   G   .   PASS    RefPanelAF=0.385864;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   255395  rs331869    C   T   .   PASS    RefPanelAF=0.980998;AN=2;AC=2;INFO=1    GT:ADS:DS:GP    1|1:1,1:2:0,0,1

Is my genotype CG or TT at this position?
Point 3

Sometimes a rsID appears twice but at different locations:

3   4942430 rs71634747  G   A   .   PASS    RefPanelAF=0.445673;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    1|0:1,0:1:0,1,0
3   4942432 rs71634747  C   G   .   PASS    RefPanelAF=0.248907;AN=2;AC=0;INFO=1    GT:ADS:DS:GP    0|0:0,0:0:1,0,0

Unless I'm really confused, this just seems wrong. If it wasn't for this weirdness, I'd seek concrete answers for #1 and #2, but it could be that they are both just a symptom of this potential error?

Oh, vr-impute@sanger.ac.uk isn't responding to my email. I have asked the above of them.

Don't know. These could be retired or inconsistent SNPs (e.g. shifted locations, or otherwise found to be different from expected). It might be that the imputation suggests that the SNP should be something other than what is observed. I can't see anything obvious from the dbSNP information that I've looked at.

a) In this specific case for rs116059446, that's a SNP that has more than two variants. One way of normalising VCF files is to rearrange them so that each line only refer to a single polymorphism from the reference sequence, so if there are two non-reference variants, the variant will get two lines in the VCF file.
b) For rs331869, this is another SNP with more than two variants recorded in dbSNP. Given that you're digging through the weeds of a huge dataset, the way of thinking about this needs to change from "Is this common" to "Is this possible". Problems are much more likely to exist with SNPs that aren't biallelic, because it indicates that something funky is going on with genetics or testing. I'll propose three possibilities for what could be happening here, in what I consider to be descending order of likelihood:

you have one non-reference allele on one of the sister chromosome copies, and another non-reference allele on another copy
you have a somatic mutation (i.e. a mutation that happened during your lifetime), and some of your cells have one variant, and others have a different variant.
the probes are tagging different regions in your genome from what is expected. Probes should be location-unique, but that's impossible to guarantee for everyone, because genetic variation is a thing that happens.

In this specific case for rs71634747, it looks like that SNP has been retired, and has merged with rs724135. It's possible that this is a multi-nucleotide variant where variants are only observed in tandem. This is an intronic variant, so seeing a 3-nucleotide variation is a little unexpected, but given that you're picking out exceptions from a very large dataset, not completely unexpected.

I think the main take-home from this is that you shouldn't ever expect results to be a perfect representation of reality. It can be useful to look at raw data to help understand what's going on, but there's too much weird stuff going on with biology, chemistry and physics for anything to be 100% accurate.
If you did exactly the same test at a later date, some of your 23andMe results will be different. If you did imputation at a later date, with a different program version, or different imputation reference dataset, your results are likely to be different.
Accuracy can be substantially improved by using different technologies to attempt to discover the same thing, and only concentrating on consistently-reported differences. You could try an Infinium Omni2.5 SNPchip, for example, or whole-genome sequencing. Unfortunately, this will lead to another problem: what should be done with the inconsistent results? The answer to that is situation dependent, and is one of the hardest challenges of bioinformatics.

Thanks for these answers. I'm not expecting data to reflect reality (I've been doing bioinformatics for too long for that ;-) I just expect data to be sane. This is just a basic sanity test, flagging things that looked wrong. I think point 1 and 3 are now tidied up (I don't think the software is behaving correctly, but I can at least understand how the software is behaving, and it's not bonkers.

I still don't understand point 2. Lets look at your potential sources for this data sanity check failure... 2b.1) makes perfect sense, but I would expect the genotype to be 0/1 and 0/1. That would be easy to interpret. That's not what I see (1/0 and 1/1). 2b.2) This may be true, but the data could never show that (imputation). 2b.3) Same, not possible in this data (imputed). Increasingly I think this call is just wrong. Not biologically wrong, just fails basic sanity checks.

Like imagine I showed you a vcf that had 3 3333 . A G . . . GT 4/4... Well, you could say, ah ha, this is a sign of a quadruplex DNA at this location in this individual... well maybe. I'd just say the VCF is broken. It's just wrong / impossible to interpret (imho).

The answer is to relax and throw away the 0.001% of stuff that doesn't make sense, lower your standards and just get on with finishing the job you set out to do without actually caring about consistency, logic, meaning or formats that are actually defined or interpretable.
Nobody cares about what these SNPs mean, how to interpret them programmatically, why there are version issues in the output of the imputation server or why it throws away ground truth in the output, leaving the user to struggle with all these issues. I mean... if nobody cares, why should I?
I've written up as detailed a description of the problems as I can here:
https://github.com/richarddurbin/pbwt/issues/52 and here: https://github.com/richarddurbin/pbwt/issues/51
and written a script to throw away all the junk here:
https://github.com/Geromics/covcheck/blob/wip/report-v2/Research/debug_imputation_results.py


text:
R help in Jupyter notebook?
I'm using this Jupyter notebook to work on R... How do I do the equivalent of ?someFun?
When I run that in the notebook, I get:
Warning message in seq_len(head.end.idx):
“first element used of 'length.out' argument”
ERROR while rich displaying an object: Error in seq_len(head.end.idx): argument must be coercible to non-negative integer

Traceback:
1. FUN(X[[i]], ...)
...

Update
When I execute docker run jupyter/r-notebook I notice that it does pop up the widget (after throwing the error). Not it isn't a popup window, but a panel at the bottom of the workbook. My notebook is 'Trusted', which I think is significant, but I don't know how to turn it off.
When I execute docker run -e JUPYTER_ENABLE_LAB=yes jupyter/r-notebook I get the error and the widget does not pop up. I'm not sure how to set or inspect the JS trust in that environment. I suspect it's related. I don't see any blocked popups in the browser.
Strange that even when it works (the widget opens at the bottom of the browser) I still get an error in the output.

If you launch like docker run -p 8888:8888 jupyter/r-notebook you will get the help page like a widget and a warning in your cell output. Check your browser, maybe is blocking pop-up elements.

Thanks for showing me what the widget looks like. I do see this now (i.e. working correctly), but only when not setting JUPYTER_ENABLE_LAB=yes

I met the same problem. R after v4.1.1 seem to have problems with repr_html.
I fixed the problem by following the quotes.

This tiny commit solves a long standing problem while trying to read help from inside jupyterlab.
In order not to wait for the PR to be merged, you may just install from @renkun-ken 's fork:
devtools::install_github("renkun-ken/repr", ref = "fix-html-filter")

Reference
https://github.com/IRkernel/repr/pull/144#event-5589475824
https://github.com/IRkernel/repr/issues/143

Thanks. I mentioned this here:
https://github.com/jupyter/docker-stacks/issues/1494

The PR is now merged and installing repr==1.1.4 resolved the issue for me.

I think error is not related to ?. This is what I see when I try ?seq:

Yes, help(seq) and ?seq are both a problem. I think it's related to trusted JS execution... I think...

source: Stack Exchange

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

sort: votes


text:
Determining a Gene's Unknown Function
I'm a beginner in this field (so pardon my rudimentary knowledge/explanation) trying to understand ways on how I could potentially identify genes of unknown function. I am trying to get as close as I can to the unknown gene's actual function.
What I've currently done, is that I've used an annotation tool to annotate my genome sample for its genes. From it, I've acquired a list of potential genes of unknown function, and I am capable of collecting its sequence. From here, I've also used BLASTn and BLASTp to align my gene sequences to find similar gene sequences.
From here, I'm absolutely clueless as to how I could proceed, since similarly sequenced genes are not always similarly functioned.
I would also like a list of other annotating tools, as well as other sequence alignment tools like BLAST, just so I could see the differences between tools.
Any help would be much appreciated.
Thank you.

Hi more detailsabout the genome, annotation and Blastp output please

Similar sequences may not always have the same function, but they're often similar enough to start with. Don't let the perfect be the enemy of the reasonably sufficient.

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


text:
about rescaled data (median)
When we have a rescaled data with median (divided by median value)(e.g.,metabolome) and perform natural log transformation, do we need to use autoscaling again to perform PCA or to compute euclidian distance?

I understand that PCA and computing Euclidian distance need scaled data. Do I need to autoscale again to these median-rescaled data (simply divided by median), Some companies give us rescaled proteome data.

I would like to know when we have median scaled data (divided by median) in metabolomics , some use natural log transformation and to compute euclidian distance and others use natural log transformation with subsequent autoscaling to compute euclidian distance. So I would like to know which methods we should use.

Whether to scale your data again in the PCA function even if you have already had some reasonable ways to scale it beforehand has no clear answer. Theoretically, there is no clear reason to do it or not. I guess you could try both and look at the resulting figure and decide.


text:
Is there a way to customize what observations are plotted in ComplexHeatmap?
So right now my observations in the dataset I've processed is stored as sample identifier columns and gene rows: an example of this would be for sample 1 and gene TRIM21 the observation is Missense Mutation;Frame Shift Deletion. I'm having issues separating the two alterations and instead of being able to show the gene in the sample having 2 different alterations, I end up having to generate a new colour classification (being that the gene in the specific sample exhibits specifically both a missense mutation and a frame shift deletion instead of being able to represent the observation with just colours and missense and frame shift deletion colours alone). Here is an a screenshot of what my matrix looks like:
As an example here is some test code generating the test plot:
library(ComplexHeatmap)

test_labels <- c("Missense_Mutation", "Splice_Site", "Frame_Shift_Del",
                 "Splice_Region","Missense_Mutation;Splice_Site", "Splice_Site;Frame_Shift_Del",
                 "Splice_Site;Missense_Mutation", "Frame_Shift_Del;Splice_Site")

plot <- Heatmap(test, row_names_side = "left",
        row_names_gp = gpar(cex = 0.7),
        column_names_gp = gpar(cex = 0.7),
        heatmap_legend_param = list(title = "Alterations",
                                    at = test_labels, labels = test_labels), show_heatmap_legend = TRUE)

plot

I've attached a screenshot of my actual plot and as you will be able to see, this ends up creating multiple different classifications:

Dear John, I am sure your post makes sense to someone, but makes little sense to me. I can see in your tags you are using Bioconductor ... that is a good start. Perhaps you can post your code and explain the data set. I see that the red is stuff you don' twant us to see

Hi Michael, thanks for responding! This is specific to a package in Bioconductor: https://jokergoo.github.io/ComplexHeatmap-reference/book/index.html. I would add my code, but I would also need to provide the exact dataset as well as the data processing I've already done, which is a decently long script. Essentially, the data set is represents samples (columns) and genes (rows), and observations correspond to what the alteration present in the specific gene in the specific patient sample. As you can see, sometimes observations are separated by ";".

This is a pretty specific question I realize, as I am using the Heatmap function in the vignette and it doesn't exactly cover how I'd be able to parse what I'd like to.

Thanks. Don't need to provide your data set ... a dumbie data set is fine - we understand patient confidentiality issues. Putting some code there would be useful.

I've updated my question with some sample code! You'll see in the legend, e.g. "Missense_Mutation;Splice_Site" is given its own category instead of being recognized separately as "Missense_Mutation" and "Splice_Site".

What color do you expect when the algorithm encounters a value that is neither Missense nor Frame-shift but Missense;Frame-shift? Are you looking for the cell to be split into 2 colors if that happens?

Because if that's the case, eyeballing will be extremely difficult with each cell looking a different height (or width) and the plot showing jarring patterns. Pick the most severe consequence and plot those alone, or use cell_fun to create a plot that IMO would be a mess and beats all reason for visualization.

What is fed into ComplexHeatmap is what it outputs! So if one of your genes has a Missense_Mutation as well as a Frame_Shift_Del and you would like to plot its value, of course you will have a category and color that is different than those of Missense_Mutation or Frame_Shift_Del, that is only natural.
What you can do is to create additional variables like your_gene_Missense_Mutation or your_gene_Frame_Shift_Del each with 0 or 1 to indicate if a particular mutation type is present or not for a particular gene (you can think of it like one-hot-encoding). With a lot of different levels, the resulting heatmap could be difficult to interpret though.
Another simpler alternative would be a using a "universal label" if your gene has multiple mutations. Basically if the cell contains a ;, its value will be changed to something like mixed_mutation_status and instead of multiple "mutation_type_1; mutation_type_2; ... like categories, you will just have one mixed_mutation_status on top of Missense_Mutation, Frame_Shift_Del, ...

The mixed_mutation_status is indeed how maftools does it - it calls it multi-hit. Although, in this case, the same mutation can be both missense and splice_site if it's a mutation at an exonic spllice site. OP should first pick among the many possible consequences of a mutation wherever possible, and then apply a logic where multiple mutations still exist in the same gene.


text:
Finding orthologs
I am interested in finding an ortholog for a particular gene in fungi.

Hi what are you goals for the research?

Yes, please [edit] your question and give us more context. What species is your original gene in? What fungi are you looking at? Are they close to your starting species? Far? What kind of orthology are you looking for? Do you require functional equivalence? Does your gene code for a protein and does this protein have any functional domains? You almost certainly want to start with blasting your protein against a translated genome (tBLASTn) but the details will depend on what you're after exactly.


text:
PCA plot did not work in single cell RNA-seq
I run plotPCA for single cell RNA-seq data, while I get error message (I use R 4.0). I attached the code and error message here. Did anyone know the reason and how to deal with it. Thanks!
library(pcaMethods)
library(SC3)
library(scater)
library(SingleCellExperiment)
library(pheatmap)
library(mclust)
set.seed(1234567)
deng <- readRDS("deng-reads.rds")
plotPCA(deng, colour_by = "cell_type2")

Error in value[[3L]](cond) : invalid subscript 'type' in 'reducedDimNames(<SingleCellExperiment>)'
'PCA' not in 'reducedDimNames(<SingleCellExperiment>)'

I also get error message when I run the code as below:
colData(deng)$tSNE_kmeans <- as.character(kmeans(deng@reducedDims$TSNE, centers = 8)$clust)

# Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), : 'data' must be of a vector type, was 'NULL'

I solve the first error for PCA, while have no idea about the second one. Does anyone have any suggestion?


text:
Variants list to aminoacids
I have the sequencing data of a Monkeypox virus sample. I have the list of variants extracted from Nextclade and I want to find the corresponding aminoacids substitution.
For example I have the variant C1092T (nucleotide annotation) and I want to automatically convert it into OPG001:D162N (aminoacid annotation).

Please could you provide a reproducible example of the problem you face?

https://www.ensembl.org/info/docs/tools/vep/index.html

Neither C1092T nor OPG001:D162N seem to follow HGVS conventions. You could try to convert them to HGVS nomenclature standards then use modules such as hgvs to get amino acid changes from nucleotide changes. You will need reference sequence identifiers though, so it might be easier to go from chromosomal level changes to AA changes, as cDNA changes depend on the transcript (which may or may not be an issue with viral transcripts).


text:
How to annotate list of proteins with KEGG orthology groups?
I've got a list of proteins, complete with Uniprot IDs (for some of them) and gene names. I'd like to group them into orthology groups based on data from KEGG, but I have no idea how to add a column with KEGG groups. I've tried poking around in the documentation for UniProt, but all I can find is this page on KEGG - no mention of how to link to KEGG data. I've also found these KEGG mapper tools (web portal and command line) but those don't seem to be what I need in order to cluster my proteins into groups of KEGG orthologues.
I've tried using R via the following:
library(KEGGREST)

kegg_ids <- KEGGREST::keggConv(
  target = "genes",
  source = paste("uniprot", uniprot_ids, sep = ':'))

But this turns up matches for only a fraction of my protein IDs... and anyway, only about half of them have those IDs - the rest are human proteins for which I only have the gene name.


text:
How can I plot multiple sequence alignment in R as a heatmap
I have many amino acids sequences in a fasta format that I did multiple sequence alignment. I was trying to plot something like binary a code as heatmap. If it had a change in the AA compared to the reference sequence it would be red, if it did not change would be yellow, for example.
I came across to msaplot from ggtree package. I also checked ggmsa package for that. But so far, I did not get what I wanted. So basically I wanted to:

change the multiple sequence alignment to a binary matrix (if the amino differs from the reference sequence, plot x, if not y)
plot as a heatmap with two colors based if the amino acid changed or not

here's some data example of what I doing, just instal ggmsa and you can try this:
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
ggmsa(protein_sequences, start = 265, end = 300, color = "Chemistry_AA",font = NULL)

Could you post your code?

I added a very simple example code present on ggmsa package

@StupidWolf answered me on stackoverflow and he's answer solved my question. I will share it here as well.
@StupidWolf:
We read in the alignment:
library(Biostrings)
library(ggmsa)
protein_sequences <- system.file("extdata", "sample.fasta", package = "ggmsa")
aln = readAAMultipleAlignment(protein_sequences)
ggmsa(protein_sequences, start = 265, end = 300)

Set the reference as the 1st sequence, some Rattus, you can also use the consensus with consensusString() :
aln = unmasked(aln)
names(aln)[1]
[1] "PH4H_Rattus_norvegicus"

ref = aln[1]

Here we iterate through the sequence and make the binary for where the sequences are the same as the reference:
bm = sapply(1:length(aln),function(i){
as.numeric(as.matrix(aln[i])==as.matrix(ref))
})

bm = t(bm)
rownames(bm) = names(aln)

The plot you see above has sequences reversed, so to visualize the same thing we reverse it, and also subset on 265 - 300:
library(pheatmap)
pheatmap(bm[nrow(bm):1,265:300],cluster_rows=FALSE,cluster_cols=FALSE)

The last row, is Rattus, the reference, and anything similar to that is read, as you can see in the alignment above, last 4 sequences are identical.


text:
Reference for the statistical methodology used in joint genotype calling
I am searching for a good reference that explains the statistical methodology used in what is variously known as Joint Genotyping, Cohort Genotyping or Multi-sample Genotyping. One with a pseudo code implementation would also be ideal.
I have found several discussing the approach in very general terms and the relative benefits of the technique but nothing with actual equations or algorithms. Failing that, if there is an open source implementation that people would consider to be a reference implementation and that has excellent code documentation that would also be great.

To the best of my knowledge, the first joint calling model was published by Le & Durbin (2011) and refined in Li (2012). They come with detailed derivation. Both samtools and older GATK are using this model. GATK now uses a different model which only works with high-coverage data but is faster to compute. There are other alternatives (e.g. in freebayes).

Here is a recent example of the derivation of the genotype model from Octopus genotype calling paper:

All samples have known ploidy and copy number, so the likelihood function of reads $R$ given genotype $g$ is:
$p(R|g) = \mathop {\prod }\limits_{n = 1}^{|R|} \frac{1}{{|g|}}\mathop {\sum }\limits_{i = 1}^{|g|} p(r_n|g_i)$
where $|g|$ is the ploidy and $|R|$ is the number of reads. The joint genotype posterior for $S$ samples is therefore given by:
$p({\bf{g}}|{{R}},{\cal{M}}_{{g}}) \propto p({\bf{g}}|{\cal{M}}_g)\mathop {\prod }\limits_{s = 1}^S p({{R}}_s|{\bf{g}}_s)$
where the genotype prior model, $\cal{M}$, is either the uniform or HWE-coalescent prior. Unfortunately, the number of genotype combinations g grows exponentially in the number of samples S, so we cannot evaluate the full posterior distribution in general. Therefore, other than for trivial cases, we first approximate the sample marginal genotype posterior distribution under the HWE model (without mutations) $(|,\cal{M})$ and use these marginal probabilities to select K genotype combinations $1,…,$ (K is user-defined) to evaluate under the full joint genotype model. The approximate posterior marginals are computed with expectation maximization.

The source code can be found here.


text:
Where can I find lncRNA expression data for different cell types?
Are there any publicly available databases providing expression data for long non-coding RNAs (lncRNAs) across cell types of multicellular organisms? Alternatively, are there lesser known UCSC tracks for this?

Example 1: I want to compare tissue-specific expression of the human lncRNA Xist across HeLa, GM12878, and HEK293 cell types [or
more].

Example 2: I want to compare tissue-specific expression of the plant lncRNA COLDAIR across cortex, phloem, and trichoblast cell
types [or more].

As long as the lncRNA is annotated you can use https://gtexportal.org/home/gene.

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

This is a fantastic resource for comparing lncRNA expression (or any annotated gene expression) across primary, non-diseased human tissues! Thanks for sharing, Mary.
Hopefully similar portals/databases will be developed in the near future for other multicellular organisms [and maybe even for cell lines?!].


text:
How to create a list of differentially expressed (DE) genes after normalization with RUVSeq?
I am using edgeR to perform differential expression (DE) analysis on a set of RNA-seq data samples (2 controls; 8 treatments). To correct for batch effects, I am using RUVSeq.
I am able to get a list of DE genes without normalization:
x <- as.factor(rep(c("Ctl","Inf"),c(2,8)))
set <- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, "DE_genelist.txt", sep=""))

Then immediately after creating the "top" object, I use RUVg to normalize:
# [...]
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg <- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, "DE_RUVg_genelist.txt", sep=""))

And I get the error:
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("SeqExpressionSet", package = "EDASeq")’ to a data.frame

I am not sure how to print the list of normalized results like I can with the unnormalized data. Ideally, I would get a file with the same format as the edgeR output (as a .csv or .txt file):
"logFC" "logCPM" "LR" "PValue" "FDR"
"COBLL1" -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
"UBE2D1" -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
"NEK7" -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
"SMC6" -2.300 5.674738981329 61.8130019860261 3.7767230643666e-15 3.4974443325016e-12

How can I get a list of genes as an output after normalization with RUVSeq?

You do the normalization before running your edgeR. The purpose of RUVg is to remove "Remove Unwanted Variation Using Control Genes". In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.
Using the example dataset in vignette:
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set <- betweenLaneNormalization(set, which="upper")

set1 <- RUVg(set, spikes, k=1)

You can look at it, it's an expression set with counts etc, not results:
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples
  element names: counts, normalizedCounts, offset
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

You run edgeR now on the results of RUVg:
design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)

This is great! Definitely answers the question. [Although the very last line topTags(lrt) will only provide a list of 10 genes, so I'll use topTags(lrt, n=nrow(set))$table for a full list].

Follow-up question for you then: What is the utility of performing the "upper-quartile" (UQ) normalization before using RUVg and then running edgeR? My thought was, "if I am using RUVg to correct for batch effects, why should I also use UQ?" Perhaps it is important?

This is what I know about UQ normalization: scaling with or without UQ before RUVg will usually lead to similar results. It is only important to include the edgeR offset -- calcNormFactors(y, method="upperquartile") -- if you decide to scale prior to RUV use. If no UQ scaling is applied, the edgeR offset should be set to 0 (zero).

The betweenLaneNormalization() in limma adjust for signal intensity difference (in arrays) or in edgeR, you try to normalize between samples such that most genes are unchanged. RUVg adjust for unseen variation on top of this, for example batch effects. So you don't want to introduce this as a kind of variation for RUVg model to pick up. Hence you do the normalization before. And as you noted, if you use edgeR,  RUV also ask you to provide the offset , which is the normalization factors.

Hope this answers your question. You can try and see for yourself what is the offset, and what RUVg does on top of it

Wonderful explanation! Thanks for your help @StupidWolf!!

I have not used this package but from your code it seems that ruvg is not a table. Instead, it is an R object, which means that you cannot use write.table. I think the results you want is stored in the object. All R objects contain "slots" of data, which can be accessed by @. If I were you, I would type ruvg@ and should be able to see which data slots are contained in the object.

This certainly gives me more direction than I had before, so I appreciate that. Just looking at all the data stored by the ruvg object (Large SeqExpressionSet (1.1 Mb)) tells me that it may take some time to find what I am looking for in there. I'll report back if I find the results I want. Thanks for your suggestion!


text:
How to subset a GRanges object based on a specific genomic window of interest?
Goal: To keep only rows of a GRanges object if their start coordinate is within 1 Mb of the start coordinate of a gene of interest.
I am aware of some functions that come close:

GenomicRanges::nearestKNeighbors() (and related distanceToNearest(), nearest(),
precede(), follow())
GenomicRanges::subsetByOverlaps()
GenomicRanges::window()

And I am aware of the UNIX command-line 'Bedtools' tools: window, closest, and intersect, however I wish to accomplish this task using R only and the above functions have not been doing exactly as I would like.
Here is an example GRanges dataset:
> gr <- GRanges(
        seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
        ranges = IRanges(start = seq(300000, 2000000, length=10), end = seq(700000, 6000000, length=10), names = head(letters, 10)),
        strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
        score = 1:10, GC = seq(1, 0, length=10), GeneSymbol = c(paste("gene", seq(1:10), sep="")))

> gr
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  a     chr1   300000-700000      - |         1  1.000000       gene1
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  e     chr1 1055555-3055555      * |         5  0.555556       gene5
  f     chr1 1244444-3644444      + |         6  0.444444       gene6
  g     chr3 1433333-4233333      + |         7  0.333333       gene7
  h     chr3 1622222-4822222      + |         8  0.222222       gene8
  i     chr3 1811111-5411111      - |         9  0.111111       gene9
  j     chr3 2000000-6000000      - |        10  0.000000      gene10
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Ideally, by specifying the name (e.g., "gene3") or start coordinate (e.g., chr2 AND 677777) of the gene of interest AND a window size (e.g., 1000000 or 1 Mb upstream AND downstream) to a function, I hope to get an output like below:
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Notice that only three rows from the original GRanges object remain in this example output because the start coordinates for all other genes were either (i) on different chromosomes or (ii) were more than 1 Mb away from the start coordinate of gene3.

I have tried various methods to isolate the rows within 1 Mb of my gene of interest with no luck:
> gr1 <- gr[grep("gene3", gr$GeneSymbol),]
> grOut <- subsetByOverlaps(gr, gr1, maxgap = 1000000, invert=T)

> grOut <- gr[mcols(gr)$GeneSymbol == "gene3" & start(gr) < 677777-1000000 & end(gr) > 1877777+1000000]

> grOut <- window(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)))
Error: 'start' must be a single number or NA

> grOut <- gr[IRanges(start=677777, end=1877777)]
Error: subscript contains out-of-bounds ranges

> grOut <- nearest(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)), select="all", ignore.strand=T)

The last attempt recorded above was likely the closest I've gotten to my desired output by storing the index of rows in the GRanges object falling within the range specified, but it presents this index in the form of a SortedByQueryHits object, not a GRanges object as I require. It seems simple enough a task, but I haven't solved it yet; any thoughts or solutions?

For those wishing to find/verify the coordinates for a gene of interest (GOI) in their GRanges object when only the gene symbol/ID/name is known: use gr[grep("<GOI>", names(gr)),] when column names have been set as gene names ... use gr[grep("<GOI>", rownames(gr)),] when row names have been set as gene names ... use gr[grep("<GOI>", mcols(gr)$<geneNameColumn>),] when a metadata column contains gene names.

Specifying the range of interest (+/- 1 Mb of start coordinate) as a new GRanges object and then subsetting the original GRanges object based on the overlap of genes with the query works:
library("GenomicRanges")

> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]
> range <- 1000000

> q <- GRanges(seqnames=chrPos,
               ranges=IRanges(start=startPos-range, end=startPos+range))
> grOut <- subsetByOverlaps(gr, q)
> grOut
GRanges object with 3 ranges and 3 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2 488888-1288888      + |         2  0.888889       gene2
  c     chr2 677777-1877777      + |         3  0.777778       gene3
  d     chr2 866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Credit to @zx8754 for inspiring this solution.

If preferred, this can be achieved with base R functions:
# First store the chromosome and start position of your gene of interest in variables
> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]

# Then add a column to the GRanges object with distances from the start position of each
# gene in the GRanges object to the start position of your gene of interest
> gr$distance <- start(ranges(gr))-startPos

# Finally subset the GRanges object based on your desired window of interest (must be
# symmetrical —— e.g., +/- 1 Mb —— for the following to work)
> grOut <- gr[seqnames(gr) %in% chrPos & (abs(gr$distance) <= 1000000)]
> grOut
GRanges object with 3 ranges and 4 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol  distance
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character> <integer>
  b     chr2 488888-1288888      + |         2  0.888889       gene2   -188889
  c     chr2 677777-1877777      + |         3  0.777778       gene3         0
  d     chr2 866666-2466666      * |         4  0.666667       gene4    188889
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

This strategy provides the added benefit of adding a 'distance-from-gene' (gr$distance) column to the dataset which can be useful for visualizations.


text:
How does DESeq2 "collapseReplicates()" function work on read counts data?
Comparing read counts from an RNA-seq experiment for two select genes before and after using DESeq2's collapseReplicates() and plotCounts() functions yields interesting results:
Before collapseReplicates() and plotCounts():

Geneid
foo1.1
foo1.2
foo2.1
foo2.2
bar1.1
bar1.2
bar2.1
bar2.2
baz1.1
baz1.2
baz2.1
baz2.2
baz3.1
baz3.2

WASH7P
6
5
0
2
1
1
8
5
0
0
0
0
0
0

SOX3
1880
1861
1950
2055
1189
1181
2415
2482
3887
3810
1851
1738
3217
3406

After collapseReplicates() and plotCounts():

Geneid
foo1
foo2
bar1
bar2
baz1
baz2
baz3

WASH7P
9.877191
2.279384
3.478891
11.613875
0.500000
0.500000
0.500000

SOX3
3189.598
3563.717
3530.486
4187.011
6473.122
7991.460
5390.721

Note: In the above comparison, there are two (2) technical replicates (1.1 + 1.2, 2.1 + 2.2, etc.) for each biological replicate (foo1, foo2, etc.) for each of three (3) conditions (foo, bar, and baz).

Comparing the tables above, it appears as though DESeq2 is NOT taking the average or sum of columns being collapsed.
It is also curious - and mildly concerning - that some very low expression (i.e., 0.5 counts) is reported for genes in the matrix of collapsed replicates when, in the original count matrix, zero (0) reads were counted as 'mapped to that gene'.
So, how does collapseReplicates() "combine counts into single columns of the count matrix" as is described in the DESeq2 vignette?

Here is the code to collapse replicates and retrieve the number of read counts for a specific gene (e.g., WASH7P) in a dataframe to be used in a count plot, getting the "after" results shown above:
dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
keep <- rowSums(counts(ddsColl)) >= 10
dds <- ddsColl[keep,]
dds$condition <- relevel(dds$condition, ref = "baz")
dds <- DESeq(dds)
countsdf <- plotCounts(dds, gene="WASH7P", intgroup="condition", returnData=TRUE)

Hmmm... alright. I have updated the question with the code (at the bottom). Let me know if you need more, like the input count matrix (cts) or annotation file (coldata).

I just realized reading the "Plot counts" section of the DESeq2 vignette that the function plotCounts "normalizes counts by the estimated size factors (or normalization factors if these were used) and adds a pseudocount of 1/2 to allow for log scale plotting", which would explain the "0.5" I observe for the genes with "0" reads mapped to those genes. This is good to know, but I now wonder if this data can be used to accurately illustrate the expression of these genes in a plot?


text:
How to remedy a DESeq2 collapsing technical replicates error?
Goal: To ensure "the sum of the counts for [my samples] is the same as the counts in the [samples] columns in ddsColl" after collapsing replicates using DESeq2.
Before performing differential expression (DE) analysis on my RNA-seq data stored in a DESeqDataSet object, I need to collapse the read counts observed for each of my technical replicates. This should be simple enough using the DESeq2 manual and vignettes available as reference.
This is how you collapse replicates:
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
# examine the colData of the collapsed data
colData(ddsColl)
# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

However, when I try this, I get the following error with no apparent resolution reported in the manual or on online forums:
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

How might I circumvent this error and continue on with the analysis?

Input files:
> head(cts)
            KO1_P1_L001.Counts KO1_P1_L002.Counts KO1_P2_L001.Counts KO1_P2_L002.Counts KO1_P3_L001.Counts KO1_P3_L002.Counts KO2_P1_L001.Counts KO2_P1_L002.Counts KO2_P2_L001.Counts KO2_P2_L002.Counts KO2_P3_L001.Counts KO2_P3_L002.Counts WT_P1_L001.Counts WT_P1_L002.Counts WT_P2_L001.Counts WT_P2_L002.Counts WT_P3_L001.Counts WT_P3_L002.Counts
DDX11L1                      1                  0                  1                  1                  2                  1                  0                  1                  2                  1                  3                  1                 7                 2                 1                 1                 0                 0
WASH7P                     124                144                128                151                102                118                 39                 41                103                 87                105                106               125               120                13                26               104               136
MIR6859-1                    8                  6                  4                  5                  4                  6                  2                  2                  7                  2                  7                  6                 6                 6                 2                 2                 2                 9
MIR1302-2HG                  1                  4                  0                  0                  0                  2                  2                  2                  2                  2                  0                  4                 1                 2                 0                 3                 2                 2
MIR1302-2                    0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0
FAM138A                      0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0

> print(coldata)
                   condition sample           run
KO1_P1_L001.Counts    longKO     l1  long_P1_L001
KO1_P1_L002.Counts    longKO     l1  long_P1_L002
KO1_P2_L001.Counts    longKO     l2  long_P2_L001
KO1_P2_L002.Counts    longKO     l2  long_P2_L002
KO1_P3_L001.Counts    longKO     l3  long_P3_L001
KO1_P3_L002.Counts    longKO     l3  long_P3_L002
KO2_P1_L001.Counts   shortKO     s1 short_P1_L001
KO2_P1_L002.Counts   shortKO     s1 short_P1_L002
KO2_P2_L001.Counts   shortKO     s2 short_P2_L001
KO2_P2_L002.Counts   shortKO     s2 short_P2_L002
KO2_P3_L001.Counts   shortKO     s3 short_P3_L001
KO2_P3_L002.Counts   shortKO     s3 short_P3_L002
WT_P1_L001.Counts         WT     w1    wt_P1_L001
WT_P1_L002.Counts         WT     w1    wt_P1_L002
WT_P2_L001.Counts         WT     w2    wt_P2_L001
WT_P2_L002.Counts         WT     w2    wt_P2_L002
WT_P3_L001.Counts         WT     w3    wt_P3_L001
WT_P3_L002.Counts         WT     w3    wt_P3_L002

My workflow:
> dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
   condition      sample           run               runsCollapsed
    <factor> <character>   <character>                 <character>
l1    longKO          l1  long_P1_L001   long_P1_L001,long_P1_L002
l2    longKO          l2  long_P2_L001   long_P2_L001,long_P2_L002
l3    longKO          l3  long_P3_L001   long_P3_L001,long_P3_L002
s1   shortKO          s1 short_P1_L001 short_P1_L001,short_P1_L002
s2   shortKO          s2 short_P2_L001 short_P2_L001,short_P2_L002
s3   shortKO          s3 short_P3_L001 short_P3_L001,short_P3_L002
w1        WT          w1    wt_P1_L001       wt_P1_L001,wt_P1_L002
w2        WT          w2    wt_P2_L001       wt_P2_L001,wt_P2_L002
w3        WT          w3    wt_P3_L001       wt_P3_L001,wt_P3_L002

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
logical(0)

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

From the output of dds and ddsColl, it seems the number of columns decreased from 18 to 9, as expected. The only difference I see between my DESeq2 workflow output is the result of matchFirstLevel (see below for expected result). Why might this be for my dataset, but not the example dataset?
Example workflow:
> dds <- makeExampleDESeqDataSet(m=12)
> dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
> dds$run <- paste0("run",1:12)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
        condition   sample         run runsCollapsed
         <factor> <factor> <character>   <character>
sample1         B  sample1       run10   run10,run11
sample2         A  sample2        run5          run5
sample3         A  sample3        run6          run6
sample4         A  sample4        run1     run1,run4
sample5         B  sample5        run7          run7
sample6         B  sample6        run9          run9
sample7         A  sample7        run3    run3,run12
sample8         B  sample8        run8          run8
sample9         A  sample9        run2          run2

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

Very fair. I wanted to include my input data for reference and reproducibility, but perhaps it makes the content look too daunting/long? I tried to clean it up a little. I think the outputs of the console commands are still necessary because they provide important info for comparison between observed and expected outcomes.

@Gawain It's best to remove lines like colnames(ddsColl) and read.table(...)  since it doesn't affect other people's understanding of the analysis, presumably it's only there for your benefit. Only include lines of code which are exactly necessary for others to understand the analysis you did.

@user438383 I see your reasoning; I have removed those lines from the content :)

sample is a character vector, so levels(dds$sample) is returning null. Presumably you mean to use dds$sample[1] here instead, which will likely solve the problem.

This worked for me! Indeed, in my dataset, sample is a character vector, and in the example dataset, sample is a factor. That would explain why the example code was not working for me. Replacing levels(dds$sample) with dds$sample[1] in the matchFirstLevel variable creation line fixed the problem.


text:
Why are Minimap2 alignments different with CIGAR generation flag?
I am using Minimap2 (v2.26-r1175) in Linux to generate a sequence alignment between the Streptomyces coelicolor A3(2) chromosome (ref.fa) and the Mycobacterium tuberculosis chromosome (query.fa). My desired output is a PAF (Pairwise mApping Format) file.
The general way to align reference and query sequences with Minimap2 is the following:
minimap2 ref.fa query.fa > approx_mapping.paf
# OR equivalently...
minimap2 -k15 -w10 ref.fa query.fa > approx_mapping.paf

You can get Minimap2 to generate custom CIGAR tags (cg:) in the PAF by adding the -c flag:
minimap2 -c ref.fa query.fa > approx_mapping_cigar.paf
# OR equivalently...
minimap2 -c -k15 -w10 ref.fa query.fa > approx_mapping_cigar.paf

Moving to an R environment [with the intention of downstream data visualization], you can quickly assess the number of observations in each object (a.k.a. rows in each dataframe) with dim():
file1 <- "approx_mapping.paf"
paf_basic <- read.table(file1, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file1, sep = "\t")))))
dim(paf_basic)
[1] 205  18

file2 <- "approx_mapping_cigar.paf"
paf_cigar <- read.table(file2, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file2, sep = "\t")))))
dim(paf_cigar)
[1] 200  24

Note the unqiue way we have to read the PAF files into R using read.table() to ensure no columns are cut off and added to new rows. This is because read.table() is designed to examine only the first five lines of an input file to determine the number of columns in the dataframe. Any columns on lines beyond the fifth line that exceed the initially determined count are appended to a new row. In the case of PAF files generated by Minimap2 with CIGAR tags, the number of columns can vary per line, which makes this necessary.
With the default settings, I get 205 observations and with the addition of the CIGAR flag, I get 200 observations.
Why do I get fewer computed alignments when the CIGAR tags are added? I would expect the alignments between the reference and query to be identical between these runs.

Comparing the PAF files (directly generated by minimap2) after sorting (sort -k3,3 -k4,4) using vimdiff reveals that the alignments computed between the default and CIGAR setting runs are largely similar, but not identical:

Note that the image above only shows the first 44 rows of both files.
Due to the small differences seen in alignment coordinates between these two dataframes (see fields $3 and $4 in the vimdiff screenshot), it is NOT trivial to determine which five observations (rows) are new alignments computed when using the CIGAR flag -c.

Plotting these dataframes confirms that the paf_basic and paf_cigar objects do, in fact, contain some different sequence alignments from each other:
library("ggplot2")
ggplot(data = paf) +
  geom_segment(aes(x = V8, xend = V9, y = V3, yend = V4), size = 2) +
  labs(title = "paf_name", x = "Reference (Mb)", y = "Query (Mb)") +
  coord_fixed(ratio = 1) +
  theme_bw()

So, it appears as though minimap2 -c produces a slightly different alignment than minimap2, which is surprising to me and begs the question: Which alignment is better or preferable?

You've shown that you have 5 extra rows when you read it into R. However, you haven't actually shown what the extra rows are. If we could see them, we might be able to guess what's going on. I can't reproduce the issue on my own. Can you rule out the trivial explanation that R's read.table is parsing the PAF incorrectly? That seems most likely given your observation about tag-only rows.

Great point! Unfortunately, there is no straightforward way to determine what the 5 extra rows are because the two data frames have no matching observations between them (checked by dplyr::semi_join(paf_basic, paf_cigar, by=c("V3","V4","V8","V9"))). Based on the graphs included in the question (a recent edit), it seems that many of the alignment coordinates are close, but not identical. Unfortunately, the dput() outputs for these data frames are too large for SO question character limits. However, I have verified that the *.paf files match the paf R objects produced by read.table().

have you looked at the two PAF files side by side in e.g. vimdiff? It should not be too hard to find the differences leading to 2 extra rows. But IDK exactly what is going on. If you've demonstrated that adding the -c option changes the alignments themselves,  that certainly sounds wrong. I'd suggest posting an issue on minimap2 github, Heng is pretty responsive. (If you do so, be sure to link to the issue here!)

Thanks for the idea Maximilian -- I hadn't heard of vimdiff before. I gave it a go and found it immensely useful. I've updated the question to include a snapshot of the two files sorted and side-by-side in the program which reveals just how similar (but not identical) the data frames are. I'll look into posting the issue on the minimap2 github. Thanks for taking time out of your day to ponder this issue and help troubleshoot with me. Cheers!

Note also that the scores are very different with and without -c. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea

Awesome analysis Colin! Interesting to see how much of an effect -c has on identity scores of different length alignments. Thanks for adding it to the conversation here.

The program's author pointed me to the Minimap2 FAQ:

1. Alignment different with option -a or -c?
Without -a, -c or --cs, minimap2 only finds approximate
mapping locations without detailed base alignment. In particular, the
start and end positions of the alignment are impricise. With one of
those options, minimap2 will perform base alignment, which is
generally more accurate but is much slower.

So, to specifically address the main questions asked here:

Why are there fewer and non-identical computed alignments when the CIGAR tags are added?

minimap2 alignment accuracy differs when running minimap2 -k15 -w10 vs. when running minimap2 -c -k15 -w10. It is actually not unexpected that the alignment coordinates in each PAF file are only similar, not identical, and an increase or decrease in the number of alignments between the outputs is simply an artefact of this accuracy discrepancy.

Which alignment (default vs. with CIGAR tags) is better or preferable?

Although both alignments report very similar sequence relationships, adding the -c flag for CIGAR tag generation is preferable to get more accurate base alignments.

Thanks for posting the answer and accepting it. Kudos to you @Gawain

Note also that the "score" can vary dramatically depending on whether -c is used or not. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea it is also alluded to briefly in the PAF.md where score is a required field but "If alignment is not available, column 10 and 11 are still required but may be highly inaccurate."


text:
Looking for genes that are expressed the same in all the samples
I have 5 samples (from a RNA seq experiment) from patients but I have no control, so since we can't look for the differentially expressed genes we want to look for the genes that have similar expression in all the sample.
We want to do that to see if some genes that can be involved are among those expressed in all the samples.
I have never done that, I do not know how to do it and I cant find any info about it, so any ideas are welcome.
Thanks!

Genes "that are expressed the same" would be "the genes that are NOT differentially expressed" so if you cannot study the latter you cannot study the former.

Please [edit] your question and give us some details. Do you already have expression data? What data do you have? What kind of format are they in? Why would you look at all genes that are similarly expressed in all samples? That won't help you identify genes related to a phenotype, not without a control.

standard deviation/mean?

I have edit the question, I hope now is more clear. thanks!

Uhm... if your patients are humans I guess you can find controls on the Internet.

If they are not humans (!!) then you don't need differential expression analysis. Just make expression analysis and then make a script to select transcripts that are equally expressed. It's only a few lines in any language. I would use python for that.

Patients are humans but It is not easy to find controls on the Internet since the protocol and the sequencing has to be done exactly the same way,

I cannot think of a way without performing differential expression analysis. Since you don't have controls, what you could do is assigning some of the samples as controls.
Let's say you have samples A, B, C, D, E and F. You can perform differential expression analysis A, B, C (i.e. treatment) versus D, E (control). If you perform this for all possible combinations (that would be 5 choose 2 which is equal to 5 choose 3 if you decide to go for 3 vs 2 splits) and filter for genes which do not show up in any DE list you end up with the genes that are expressed similarly in al of these samples.


text:
Performing PCA for the samples and for the genes
I have 10 samples from a RNAseq experiment (5 control, 5 disease), I have performed a cluster analysis for the samples and for the genes (4000 genes aprox) to see how they cluster (to see which samples are similar and which genes have similar expressions).
I was wondering if in terms of statistics it makes any sense to perform a PCA of the samples instead of the genes?
thank you in advance!

I think doing PCA on RNAseq samples is something people do - you can probably just use the prcomp function in R. However, I'm not sure how much information you can get in a PCA with just 10 samples.

People usually do PCA on RNAseq samples, but they use genes as variables. But I want to do it using samples as variables, and I am not sure if that makes sense in terms of statistics.

Why do you specifically want to do a PCA using samples? Can you do it on genes and plot the loading of samples? Since you have thousands of genes and only 10 samples, it doesn't seem you can do it on samples.

First of all, I am not a statistician but have been performing RNA-seq analysis for a while, here is my take:

I think performing PCA on the samples makes sense in terms of mathematics. What I mean is you can compute eigenvectors and subsequently principal components either on your genes or on your samples.

I don't think performing PCA on the samples makes sense in terms of statistics in RNA-seq context (and probably in any other context): Your genes (features) describe your samples and not the other way around. So the computed PCs should describe/explain your samples and not the other way around. You would not perform PCA on your samples (individual cells in scRNA-seq) even in the case of scRNA-seq, where number of samples sometimes exceeds the number of genes, there PCA is performed on genes and the resulting few PCs are used for downstream analysis, usually clustering and visualization, both used to infer about samples (cells).

If you use, say, R's prcomp, you can use my_pca_object$x to get PC coordinates for the samples, and my_pca_object$rotation to get how much each gene contributes to each PCA.  (You might need to transpose your transformed count file to get it to work right)

It is quite possible to make a PCA plot with the samples. Follow this tutorial : https://hbctraining.github.io/DGE_workshop/lessons/03_DGE_QC_analysis.html


text:
Pathway enrichment analysis with or without custom background
I am doing a pathway enrichment analysis of a sets of proteins.
Specific cell types (cardiospheres) were used to perform the experiment and for the quantification of proteins was used mass spectrometry.
I am using gProfiler to perform pathway analysis, should I use custom background of all of those proteins identified from the experiment? or should I used without custom background?
Thank you!
UPDATE: Also, I am doing an enrichment analysis of the interactome of the set of proteins of interest. In that case should I used a custom background for the enrichment analysis or not?
Thanks!!


text:
Easy way to create gene plot in gviz?
I would like to create a simple plot of a given human gene containing the main features, e.g. exons, introns, ORFs, transcripts etc. I understand that gviz  is a popular choice for this kind of thing, but most tutorials are massive overkill for what I want to do (I don't need fancy features such as importing BAM files and plotting coverage).
I was wondering whether someone can point me towards some sample code for creating a basic gene plot. It probably would be good to start by creating a dataframe containing all the features of my gene - is there a handy package that can do that?

If Python is an option, you could check BioPython and GenomDiagram. This should get the job done.

I assume you are trying to make something like a "gene arrow map" and [gggenes][1] can help you (with some little tweaks) in case you are familiar with ggplot2.


text:
How do I quantify a specific somatic variant?
I'm working with targeted Illumina sequencing data generated with DNA from diseased and healthy tissue (this is an age-related disease and is not cancer/neoplastic). My hypothesis is that the diseased tissue might contain specific somatic mutations in one particular gene (let's call it GENE), while the healthy tissue will not. The reasoning behind this: The same mutations cause autosomal-dominant versions of my disease of interest when present in the germline (my diseased tissue is from sporadic disease). My coverage varies, but is in the range of 1000x - 50'000x. My sequencing protocol employed UMIs.
I have run Mutect2 in default settings and with lower thresholds (--tumor-lod-to-emit) and haven't found any calls that correspond to these pathogenic mutations, so my results are negative. Nevertheless, I have reason to believe that my pathogenic mutations might be present at very low VAF (around 0.05-0.1 %), based on other experiments.
My question: How do I quantify the number of reads that correspond to a particular variant of interest? Let's say that base 2000 in GENE is usually an A, whereas the mutation A>T would be pathogenic. How do I count the number of Ts in position 2000?
As I'm somewhat new to bioinformatics, I'm also very open to any other ideas. Thanks!

Edit: @terdon's approach worked perfectly in cases where I used a VCF for --alleles that had been generated as an output of Mutect2, such as the following one:
VCF produced by Mutect:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Curiously, if I used a manually produced file such as the following, no output was produced, unless the variant of interest was heterozygous:
Manually produced VCF:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 rs1799990   A   G   .   .   .

Output using the manual VCF, heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  heterogyzous_sample
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Output using the manual VCF, non-heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  novariant_sample

Manual inspection of this non-heterozygous sample revealed that it contained 1 read with the ALT variant, so this doesn't seem to depend on the variant being absent.

I have a little function that can report this by using samtools to look at the read pileup. Um. It isn't the most elegant piece of code, but it does sorta do the job:
variantSupport ()
{
  genome=$1
  bam=$2;
  range=$3;
  samtools mpileup -f "$genome" -r "$range" "$bam" | cut -f 5 |
    tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c |
    awk 'BEGIN{tot=1}$2~/,/{ $2="." }
         {
           tot+=$1;
           c[$2]+=$1
         }
         $2~/[a-zA-Z]/{ varReads+=$1; }
         END{
           for(i in c){
             print i,c[i]
           }
           printf "%s of %s reads (%.1f%%), support the variant.\n",varReads,tot,varReads*100/tot
   }'
}

Add the function to your ~/.bashrc or just paste it directly into a terminal and you can then do:
$ variantSupport /path/to/genome/hg19.fa foo.bam chr13:28607916-28607916
[mpileup] 1 samples in 1 input files
A 915
* 1
. 93
915 of 1010 reads (90.6%), support the variant.

This is telling me that at position 28607916 of chr13, I have 1 read supporting a deletion (*, see https://www.htslib.org/doc/samtools-mpileup.html for an explanation of the symbols), 93 reads supporting the reference residue (T, in this case) and 915 supporting an A, so a T=> A variant.
However, a better approach is to use a variant caller and force it to call the variant you want. You can then investigate the support for it. With mutect2, you can use the --alleles option to force the calling and then use the -L option to limit calling to only the specific region. If you give the same vcf file for both options, that will return results for the variants you want and only those. For example, given this file foo.vcf.gz (which needs to be compressed using bgzip and indexed using tabix):
$ zcat 1.vcf.gz
##fileformat=VCFv4.3
##reference=/genomes/hg19/hg19.fa
##FORMAT=<ID=AB,Number=A,Type=Float,Description="Allelic Balance">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Net Genotype quality across all datasets, calculated from GQ scores of callsets supporting the consensus GT, using only one callset from each dataset">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr13,assembly=hg19,length=115169878>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG003
chr13   28607917    .   G   C   .   .   .   GT:AB:DP:GQ 1/1:0.6:395:318

You can then do:
gatk Mutect2 \
  -R /genomes/hg19/hg19.fa \
  --alleles 1.vcf.gz  \
  -L 1.vcf.gz \
  -I foo.bam \
  -O out.vcf

This will produce an out.vcf like this (I am excluding most headers for clarity):

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  331970
chr13   28607917    .   G   C   .   .   DP=1107;ECNT=1;MBQ=30,0;MFRL=353,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-2.716e+00 GT:AD:AF:DP:F1R2:F2R1:SB    0/1:1043,0:9.616e-04:1043:518,0:496,0:989,54,0,0

I can see this is not a real variant since the value for AD, the number of reads supporting the reference and the variant is 1043,0, meaning that 1043 reads support the reference and none support the variant.
This way, you can get results for your variants of interest and see if there is any evidence supporting them.

Thanks! Your answer does what I needed. However, I have noticed that the Allele Depth numbers in the VCF (out.vcf, in your code) are lower than expected. For instance, I got REF=105 and ALT=104, whereas the numbers are  REF = 436 and ALT = 380 if I extract this information with bcftools and set -d to a high number. I need the total read numbers, so how can I edit your command to achieve this?

@Nereus I don't think you need the total, actually. The ones that are not counted by my function are the the ones  samtools skips by default (see https://www.htslib.org/doc/samtools-mpileup.html; for example, it won't count bases with a fastq quality under 13), and the ones not included in the AD value are those the variant caller filtered out. You can't just use all reads, you must only take into account the ones that pass some basic QC filtering otherwise, you just get nonsense. If you see such huge differences, I would investigate the quality (both mapping and fastq) in those regions.

Sounds good, thank you for your advice

may I ask how you produced the 1.vcf.gz file in your example?

@Nereus from one I had. I just looked at a bam file I had lying around and made the vcf point to a variant in that file and ran the command I showed you.

I see. Would it also be possible to create a custom VCF if I haven't got one available, e.g. with data from Ensembl? Thank you for your time!

@Nereus um, sure, why not? I mean, a vcf is just a text file, you can make one manually if you want.

I have noticed that your approach produces the desired result if the variant is actually present; however, when the result is not present an empty VCF file is created (I'm adding an Edit above)


text:
Error pasing the following line in pdb
I converted an .sdf file containing around 50,000 anti-viral agents to individual .pdb files using OpenBabel. However, when I tried opening the .pdb files with Autodock I kept getting the error below. What do you think is causing the error.
Thank you for your time and patience.
Error Message:

“ Error pasing the following line in pdb:
HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."

You have a firewall issue, it is needing to connect outside your computer.

I have added autodock.exe and autogrid.exe as exceptions to my firewall and yet I am still facing the same issue.
Here’s is exact the error message,
“ Error pasing the following line in pdb:
HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000 ."
What do you think I should do to overcome the error. Thank you.

Well at least it is now an incomprehensible error - but it still doesn't work :-(  .. some of the pdb gurus are needed here, its a very ancient format BTW.

A PDB file has strict spacing rules as documented here.
The coordinates should be:
31 - 38       Real(8.3)      x             Orthogonal coordinates for X.
39 - 46       Real(8.3)      y             Orthogonal coordinates for Y.
47 - 54       Real(8.3)      z             Orthogonal coordinates for Z.

So trying this on your line:
'HETATM 1 C UNL 1 52355.00975659.051 0.000 1.00 0.000'[30:39]

I get 9.051 0., which is not a number.
The whole line is formatted incorrectly however.
So as you say you are using Autodock why not pass a mol2 to prepare_ligand4.py instead which is an acceptable format?
Alternatively, you could convert your multi-entry sdf to pdb files via RDKit in Python:
suppl = Chem.SDMolSupplier(sdfile)
for mol in suppl:
      # mol = Chem.AddHs(mol) # has hydrogens?
      # AllChem.EmbedMolecule(mol) # is 3D?
      AllChem.ComputeGasteigerCharges(m) # sdfs do not have partial charges
      Chem.MolToPDBFile(m, f"{mol.GetProp('_Name')}.pdb")

Hello Mr. Ferla,
Thank you for getting back to me.
Is there any option in Open Babel that will let me make the required changes so that the spacing rules aren't violated.
Best Regards,
Shahbaaz.

@Shahbaaz , amended.
Parenthetically no need for formalities. Also @ plus username is generally the preferred way to address folk here, you might get titles  wrong ;) .

Also, this other Q about Autodock protocol may be useful: https://bioinformatics.stackexchange.com/questions/11471/protein-ligand-docking-how-to-convert-protein-pdb-to-protein-maps-fld/11479#11479

Thanks a ton for the help

Good morning @Matteo Ferla
Can you briefly summarize how I can use prepare_ligand4.py or RDKit. I have looked for tutorials online but couldn't come across any. Thank you.

First, prepare_ligand4.py is for Autodock 4, not VINA. The two are different and fill different roles: complicated but adaptable vs. easy but blackbox. I am assuming you are using the former as you did not specify VINA. Here is again the link to documentation.

I need help with the following line,
"pythonsh [yourpath or .]/prepare_ligand4.py -l ligand.pdb [options]"
considering the path for mg tools is the following, "C:\Program Files (x86)\MGLTools-1.5.6" how would I alter the above line including the mg tools file path in it

If you use finder (whatever it is called on a Windows) and explore the stuff within C:\Program Files (x86)\MGLTools-1.5.6 you'll find the script in somewhere like autodock_tools/utilities

I have tried everything to split the multi entry sdf file to pdb files but have failed. Could you split the sdf file to individual files for me as I have been unable to do so despite multiple attempts at it.
Thanking you in anticipation,
Dr. Hussain.

If it is just for training, why not redock ligands from already bound structures one at the time. There are lots of inhibitors to the COVID protease in the PDB and the ligand can be downloaded as sdf, pdb or mol2. Make an apo, place the ligand (don't keep the original there) and then you can see what RSMD you get versus the actual solved one.
Once there you can try your hand at virtual library screening.


text:
How do I set a neural network to loop multiple times and average the resulting values?
I have a script in R/RStudio which creates random datasets of binomial variables, feeds them through a neural network, and calculates their likelihood ratio statistic and deviance. I'd like the script to loop with the seed incremented by 1, and collect the value from each run into a data frame, from which the values get averaged. For example, if the script returns 5, 10, and 15 on the first, second, and third iterations of the loop, I'd like the script to create a data frame with those values and then calculate their mean.
Here is my script, annotated:
library(neuralnet)

x=1021

#change nrow to number of times code is looped
result <- data.frame(matrix(nrow = 1000, ncol = 2))
colnames(result) <- c("LogL", "D")

repeat{

set.seed(x)

#2 inputs, max number of terms 3
x1=rbinom(2000,1,0.5)
x2=rbinom(2000,1,0.5)
y=rbinom(2000,1,0.5)

data=data.frame(x1,x2,y)

head(data)

#hidden units = 2
nnet=neuralnet(y~x1+x2,data,hidden=2,err.fct = "ce",linear.output = F,likelihood = T)

predictions=compute(nnet,data[,1:2])

predictions$net.result

#likelihood ratio stat
argument=y*log(predictions$net.result)+(1-y)*log(1-predictions$net.result)

head(argument)

#sum of likelihood ratio stat
LogL=sum(argument)

LogL

D=-2*LogL

D

result[x-1022, 1] <- LogL
result[x-1022, 2] <- D

result

#results <- data.frame(LogL, D)

x=x+1
if (x==2021) {
  break
}

print(x)

}

result

My repeat loop works, but a thousand iterations, for example, takes multiple minutes. Is there a more efficient way to code such a loop?

Your repeat termination conditions will never be satisfied. Every round, x will be set to 1021, incremented by 1 and compared to 1023, which will never be true. How will the repeat ever stop?

There are other approaches to parameterisation in ANN. @RamRS question on terminating recursion is definitely a good question.

Thanks, edited that and it works now, but it's super slow. I'll update the question!

The concern with your approach is called 'leakage' because you are parameterising the same data set that you are training. This can easily lead to overtraining is in ANN (artificial neural networks) and is not cool. Overtraining is where the ANN very tightly wraps around the training data set that it performs poorly on any other data set. In other words it is overtightened.
The way they get around it is to split the data set usually 60:20:20 into a training split, parameterisation (validation) split and test (testing) split. If the parameterisation split gets mixed with the training set - leakage happens and overtraining results. Normally the data is trained, and the trained model parameters are assessed/optimised in the parameter split. The trained data and optimised parameters are then tested in the test split.
As long as you are aware of the approach you are using, that I guess that is ok. If you separate your training and parameterisation (validation) data that would be cool.


text:
How to design DESeq2 LRT model with individuals nested in 2 levels?
We have a complicated experimental design that we would like to perform LRT analysis for. Our main goal is to discover significant genes for the "Injection:Social" interaction term across the entire dataset by removing it from the LRT reduced model, and as a bonus we are also interested in discovering significant genes for that interaction term for each respective brain region.
Sample  Injection   Social  Region  Individual  ind.n
HY06    L   ISO HY  S06 S1
NST6    L   ISO NS  S06 S1
TN06    L   ISO TN  S06 S1
HY08    L   ISO HY  S08 S2
NST8    L   ISO NS  S08 S2
TN08    L   ISO TN  S08 S2
HY30    L   KF  HY  S30 S1
NST30   L   KF  NS  S30 S1
TN30    L   KF  TN  S30 S1
HY32    L   KF  HY  S32 S2
NST32   L   KF  NS  S32 S2
TN32    L   KF  TN  S32 S2
HY64    L   KFC HY  S64 S1
NST64   L   KFC NS  S64 S1
TN64    L   KFC TN  S64 S1
HY65    L   KFC HY  S65 S2
NST65   L   KFC NS  S65 S2
TN65    L   KFC TN  S65 S2
HY19    L   NF  HY  S19 S1
NST19   L   NF  NS  S19 S1
TN19    L   NF  TN  S19 S1
HY24    L   NF  HY  S24 S2
NST24   L   NF  NS  S24 S2
TN24    L   NF  TN  S24 S2
HY05    S   ISO HY  S05 S1
NST5    S   ISO NS  S05 S1
TN05    S   ISO TN  S05 S1
HY12    S   ISO HY  S12 S2
NST12   S   ISO NS  S12 S2
TN12    S   ISO TN  S12 S2
HY31    S   KF  HY  S31 S1
NST31   S   KF  NS  S31 S1
TN31    S   KF  TN  S31 S1
HY34    S   KF  HY  S34 S2
NST34   S   KF  NS  S34 S2
TN34    S   KF  TN  S34 S2
HY62    S   KFC HY  S62 S1
NST62   S   KFC NS  S62 S1
TN62    S   KFC TN  S62 S1
HY63    S   KFC HY  S63 S2
NST63   S   KFC NS  S63 S2
TN63    S   KFC TN  S63 S2
HY04    S   NF  HY  S04 S1
NST4    S   NF  NS  S04 S1
TN04    S   NF  TN  S04 S1
HY20    S   NF  HY  S20 S2
NST20   S   NF  NS  S20 S2
TN20    S   NF  TN  S20 S2

My first attempt was building simple full (m1) and reduced (m2) models that gets directly at our question of interest but doesn't control for nested individuals.
m1 <- model.matrix(~ Region + Social * Injection, colData_filt)
m2 <- model.matrix(~ Region + Social + Injection, colData_filt)

We want to control for individual/batch effects, which is nested within both "Injection" and "Social" but not region, as we have three brain regions per individual. I followed the example in the DESeq2 manual for creating a term (ind.n) distiguishing individuals nested within groups, but now I'm not sure how to create the full and reduced model given that I have one more level than the example.
I've tried a really elaborate full model (m1) with the interaction term of interest (Injection:Social) removed for the reduced model (m2), but I'm not sure this is correct based on our design.
m1 <- model.matrix(~ Injection + Injection:ind.n + Injection:Social + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)
m2 <- model.matrix(~ Injection + Injection:ind.n + Injection:Region + Social + Social:ind.n + Social:Region + Region, colData_filt)

I'm assuming this is wrong, but even if this was by some miracle the correct formulation, would there be a way to extract genes that explain the "Injection:Social" interaction term for separate brain regions?
As a work-around, I subsetted the data by region and ran three separate LRT analyses for each subset and compared the results. While this simplified the model to look like the first example above, I worry that we lose some power by ignoring the fact we have multiple brain region samples from single individuals across the dataset.
Any guidance is much appreciated. Thanks in advance

From what I can gather, you want to account for the effect of individual, nested within region. That is, you want to see after accounting for these, is there a consistent effect for Injection:Social across all conditions.
So you set up the model like this:
m1 <- model.matrix(~ ind.n*Region + Injection + Social + Injection:Social,data=..)

The last term should be Injection:Region and you can just use the waldTest (default) in DESeq2 for this term.
What does the terms do? ind.n*Region is the equivalent of ind.n + Region + ind.n:Region , and with this you effectively get an effect for every region in every individual.
Why don't we need the Injection:ind.n or Social:ind.n or Social:Region. These terms indicate the effect of Injection or Social can vary by individuals or regions. Most likely introducing too many parameters when you are interested in a common effect. Also you do not have the replicates or samples to distinguish this effect from region or other effects.
Since you provided an example we can run DESeq2 and you can see how the results look like:
mat = counts(makeExampleDESeqDataSet(n=1000,m=48))
dds = DESeqDataSetFromMatrix(mat,df,~ ind.n*Region + Injection + Social + Injection:Social)

dds = DESeq(dds)

resultsNames(dds)
 [1] "Intercept"            "ind.n_S2_vs_S1"       "Region_NS_vs_HY"
 [4] "Region_TN_vs_HY"      "Injection_S_vs_L"     "Social_KF_vs_ISO"
 [7] "Social_KFC_vs_ISO"    "Social_NF_vs_ISO"     "ind.nS2.RegionNS"
[10] "ind.nS2.RegionTN"     "InjectionS.SocialKF"  "InjectionS.SocialKFC"
[13] "InjectionS.SocialNF"

The terms you need are "InjectionS.SocialKF","InjectionS.SocialKFC", "InjectionS.SocialNF", and you can look at each of them:
head(results(dds,name="InjectionS.SocialNF"))
log2 fold change (MLE): InjectionS.SocialNF
Wald test p-value: InjectionS.SocialNF
DataFrame with 6 rows and 6 columns
              baseMean     log2FoldChange             lfcSE               stat
             <numeric>          <numeric>         <numeric>          <numeric>
gene1  9.9811166787259   1.25304112986447 0.819806376919295   1.52845984752303
gene2 30.3449455820337 0.0329442893152027 0.705199688255367 0.0467162562092241
gene3 3.83223545055379    1.0281136369045  1.64095596190233  0.626533350543196
gene4  11.232305747171  0.595738624408923   0.8243883031544  0.722643227868976
gene5 6.70950627004097  0.756449993378065   1.0631622863378  0.711509430967263
gene6 26.1431134888287 -0.854784518963918 0.625714541243558  -1.36609342219393
                 pvalue              padj
              <numeric>         <numeric>
gene1 0.126398405431826 0.978671002658464
gene2 0.962739373909937 0.999897888026606
gene3 0.530965168963838 0.978671002658464
gene4 0.469899103018657 0.978671002658464
gene5 0.476768608734069 0.978671002658464
gene6 0.171909642630577 0.978671002658464

As mentioned, you can do a LRT if you want to test all the Injection:Social term interactions terms at one go, that is, the null hypothesis is that all of them are zero:
dds = nbinomLRT(dds,reduced=~ ind.n*Region + Injection + Social)
results(dds)

Usually the individual terms are more intuitively than testing all of them are zero, but you might have a special need for this.

Thank you, I will try this soon. Just to clarify, an LRT would still be appropriate for the larger question of finding the effect of the "Injection:Social" term and you suggest Wald for determining the region-specific impacts? Would it also be valid include a term like "Injection:Social:Region" in the suggested Wald model formulation in an attempt to explore the interactive effects of "Injection:Social" by region, or would that require more 3-way interaction terms with "ind.n"?

@jfaberha, thanks for providing the data, i have updated my answer, hopefully its clearer now. You can include "Injection:Social:Region" only if you include Region:Social. My point is this, do you think there is such an effect?

I would visualize this with a pca or correlation matrix to see where the effects are coming from before applying them onto the model. No need to have a super complicated model...

To answer your question regarding the "Injection:Social:Region" term, I will do some more research to see if this has an appreciable effect in the dataset. One question posed by the PI of this study is - what are the genes responsive to the "Injection:Social" interaction term in each of the different brain regions and are they different? I was guessing that maybe a 3-way interaction term could address this in the binomial Wald analysis, but I agree that adds excessive complexity to the model. That's said, I'm not sure how to get region-specific interactions without running separate analyses.

hmmm what is the main question? it sounds like everything under the sun. So it is indeed Injection:Social with different brain regions, you are better off running separate analysis for each region. It is easier for you to explain and honestly you don't gain power by analyzing all together

The main question is regarding the Injection:Social term across the dataset, with the brain region-specific interactions being less important follow-up. That said, I will move forward with LRT for the primary question and separate analyses for the region-specific analyses. Thank you for all the help!

sounds like a good plan!

Hi @jfaberha remember to upvote SW, LRTs are powerful statistics but never trivial to implement and that is a good answer.

in LRT can we use reference level which I use in pairwise comparison i have  a made a similar query Since i m not sure  if I use LRT which would be my reference in a reduced model


text:
Multiple input & output files?
I'm trying to print out the measured distances for each .pdbqt file, and then write a new output file for each .pdbqt, using the original name of the .pdbqt file as the new .txt file name. Below is the script.
I'm using python3 & I'm trying to do this via Linux CLI.
from pymol import stored
stored.a, stored.b, stored.c, stored.d = [],[],[],[]
cmd.iterate_state(1, "r55_01///UNL`1/N*", "stored.a.append(ID)")
cmd.iterate_state(1, "r55_01///UNL`1/O*", "stored.b.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASP`113/OD2", "stored.c.append(ID)")
cmd.iterate_state(1, "r55_01//R/ASN`312/ND2", "stored.d.append(ID)")

f=open("tt.txt", 'w')

for a in stored.a:
  for c in stored.c:
      d1=cmd.distance("id %s" % a, "id %s" % c,cutoff=5)
      f.write("%s D113 %s\n" % (a,d1))
      print(d1)

for b in stored.b:
  for d in stored.d:
      d2 = cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
      f.write("%s N312 %s\n" % (b,d2))
      print(d2)

f.close()

My folder looks something like this
.
├── dst6.txt
├── sp
│   ├── r55_01.pdbqt
│   ├── r55_02.pdbqt
│   ├── r55_03.pdbqt
│   ├── r55_04.pdbqt
│   ├── r55_05.pdbqt
│   ├── r55_06.pdbqt
│   ├── r55_07.pdbqt
│   ├── r55_08.pdbqt
│   ├── r55_09.pdbqt
│   ├── r55_10.pdbqt
│   ├── r56_01.pdbqt
│   ├── r56_02.pdbqt
│   ├── r56_03.pdbqt
│   ├── r56_04.pdbqt
│   ├── r56_05.pdbqt
│   ├── r56_06.pdbqt
│   ├── r56_07.pdbqt
│   ├── r56_08.pdbqt
│   ├── r56_09.pdbqt
│   └── r56_10.pdbqt
└── t6.py

THANKS IN ADVANCE!!

You made the hardest work, you were not far to achieve what you wanted ;)
This script.py should work after you made it executable (in a terminal type: chmod u+x script.py), just type ./script.py in your terminal in your working directory (so here you should be in your "sp" directory). Note that I got inspired by: https://www.biostars.org/p/42474/ to know how to launch pymol within the script and also that I would personally not use pymol for doing this kind of analysis/data extraction from a pdb/pdbqt file but rather use a dedicated analysis package (there are plenty like MDanalysis). I personally use pymol (pymol script) for visualization purpose.
#!/usr/bin/env python

import __main__
__main__.pymol_argv = ['pymol', '-qc']  # Quiet and no GUI

from pymol import stored
import pymol, os

pymol.finish_launching()

#put what you've done into a function
def myfunc(file):
    spath = os.path.abspath(file)
    sname = spath.split('/')[-1].split('.')[0]

    pymol.cmd.load(spath, sname)

    stored.a, stored.b, stored.c, stored.d = [], [], [], []
    pymol.cmd.iterate_state(1, "%s///UNL`1/N*" % sname, "stored.a.append(ID)")
    pymol.cmd.iterate_state(1, "%s///UNL`1/O*" % sname, "stored.b.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASP`113/OD2" % sname, "stored.c.append(ID)")
    pymol.cmd.iterate_state(1, "%s//R/ASN`312/ND2" % sname, "stored.d.append(ID)")

    f = open('%s.txt' % sname, 'w')

    for a in stored.a:
        for c in stored.c:
            d1 = pymol.cmd.distance("id %s" % a, "id %s" % c, cutoff=5)
            f.write("%s D113 %s\n" % (a, d1))
            print(d1)

    for b in stored.b:
        for d in stored.d:
            d2 = pymol.cmd.distance("id %s" % b, "id %s" % d, cutoff=5)
            f.write("%s N312 %s\n" % (b, d2))
            print(d2)

    f.close()

mypdbqt = os.popen('ls *.pdbqt').read().split()

for i in mypdbqt:
    myfunc(i)

pymol.cmd.quit()