text:
HTSeq-count is returning 0 for every gene, instead of expression value
(Ported from stack overflow)
I'm trying to summarize gene count using htseq-count; and it's returning 0 counts at every gene. I'm not sure what I'm doing wrong; I think it has to do with the gene flag I'm using.
I've tried using the GTF for Arabidopsis TAIR 10: which I got from: ftp://ftp.ensemblgenomes.org/pub/release-46/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.46.gtf.gz as well as the GFF file available on Araport.
I'm using various BAM files I know to have gene expression values.
I've tried both locally and on galaxy
The local commands I've tried so far are:
python3 -m HTSeq.scripts.count --idattr=gene_id #bam #gtf > summary.txt
python3 -m HTSeq.scripts.count -i gene_id #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i locus_tag #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i ID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i SeqID #bam #gff > summary.txt
python3 -m HTSeq.scripts.count -i Parent #bam #gff > summary.txt

I've done the same on galaxy; inputing the -i value in.
Any help as to where I'm going wrong would be great.
Thanks!

hey @sdshinghal, check a few things. 1. ensure bam file is ok (you can look at it), 2. like the answer below, check your gtf matches the reference genome

Maybe the contig names in your BAM files don't match those in the GFF file? Ensembl uses contig names like 1, some other reference sources use chr1.
Generally, for ensembl GFF files, HTSeq's default value for idattr should be fine though.


text:
polishing assembled genome to QV50 value
I am trying to assemble some genomes from pacbio reads using HGAP3 followed by quiver polishing. However, even after multiple rounds of quiver polishing, I cannot attain quality value 50 (QV50) for a genome.
QV50 seems to be recommended and I am able to get upto QV30 only. I have done polishing for almost 15 times already. I have tried polishing with different parameter changes like subread length cutoff and blasr maximum divergence.
Can anyone suggest me on how I can achieve an assembled genome of QV50 value?
Thank you

15 times is too many times. If you were going to get quiver to get you good QVs it would have happened after 2-3 rounds probably.
Some questions:

Why are you using quiver? Are your reads from RSII chemistry? Arrow is recommended for Sequel reads.
What is your dataset coverage? If you do not have very high coverage (I think >30X is usually a minimum for a eukaryotic genome) it may be that you are just not going to get exactly the values you want because assembling low coverage data is hard.
Have you tried other polishing tools? They might be able to get you better value. I am not an expert but some people like to use Racon or Apollo.
While I'm not an expert, looking into how the QVs are calculated can be helpful. I looked at the quiver docs and it looks like if the reads map ambiguously then QV is low. So if you for example have lots of haplotigs, then mapping might be ambiguous for a lot of reads, and therefore this will depress your QV values. This is therefore a problem that can come from the starting assembly and not just the polishing. You may want to go back and take another look at your starting assembly.

source: Stack Exchange

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

sort: votes


text:
Optimum parameters for genome assembly via canu
I have been trying to assemble my genome via canu. Canu gave me about 200 contigs using default parameters. My genomesize is ~9Mbp. Can you recommend some parameters I can tweak to decrease the contigs number?

You can also try other de novo assembly tools like SMARTdenovo to see if they output results more in line with what you expect

The canu parameter reference guide is always a good place to start.
It would help to know a little bit more about what your assembly looks like. E.g.

what is your n50, n90...?
What is your dataset coverage?
Have you run any tools like quast to evaluate quality?
What is your input data type (probably PacBio)? Read length distribution?
What is your expected chromosome/replicon number?
Have you compared the assembly to any other data, e.g. genomes of close relatives? (quast can do this)

Basically, why do you think that 200 contigs is too many contigs? Sometimes you just don't have enough coverage to get more contiguous/"nice" assemblies. That's not to say it's not possible, of course, just that we don't have enough information to evaluate.
It's possible that some of those contigs are small and artefactual, so you might get rid of them with just length/coverage filtering. You can do some of this in canu itself (see "Output Filtering" section).

It is pacbio reads. It results in a contig which is circular and ~9.5 Mbp. Other contigs are either 2Mbp or lesser. They do not seem to be plasmids as well from annotations, and their coverage is also good, ~20X (for the 9.5Mbp contig, its ~30X). Also, other contigs are not repeats of the 9.5 Mbp contig. The longest contig is similar to the reference genome from mauve alignment and dot plot.
I am confused if even the 2Mbp contigs could be artefactual?

@menten well, to me it sounds like your circular contig is your reference genome, congratulations. the fact that you get other contigs in addition to it suggests to me 1) you have contamination. 2) the other contigs are haplotigs (contigs from alternate haplotypes or clones in your culture) of your genome. I would suggest doing an all-by-all alignment of your contigs via minimap2, which should identify any haplotigs. do you know the ploidy? if it's diploid and not homozygous, well, then you likely have haplotigs with high-coverage pacbio.

(i am guessing that you are working on bacteria in which case I know high ploidy may be an odd question, but you never know)

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


text:
Running MaxQuant on Linux
I'm trying to run MaxQuant on a linux laptop, hower I'm constantly running into problems with MaxQuant crashing at different stages
I tried to use

mono 6.8
mono 5.4.1

on

centos
ubuntu

each with MaxQuant

1.6.3
1.6.7
1.6.11

(everything using docker), but I had no success so far.
I also tried to use the ahead of time compilation of mono, also without sucess.
MQ mostly stalls at the "preparing searches" stage.
Has anyone managed to get MQ running properly on linux and can tell me, which linux distro/version, mono version and MQ version you are using?
Any other help is of course also welcome

There is a conda package for maxquant that should work on most platforms (I've used it on CentOS and I think it's been used on Ubuntu as well). There's also a docker container version of that that we (Bioconda) have posted to the biocontainers project that should also work.
Please note that large parts of a maxquant analysis are single threaded, regardless of what you specify. We've also had single analysis runs lasting over a month (even when using 90 threads), so the program can be very slow if your search space is large.

I have been running MaxQuant quite frequently on Linux. I would use conda to create a dedicated environment and install mono.
# create the environment
conda create -n maxquant -c conda-forge mono
# activate the environment
conda activate maxquant
# run any maxquant version
mono /path/to/maxuant/MaxQuantCmd.exe mqpar.xml

This installs a conda environment with mono which is a .NET framework for Linux. I had problems in the past running certain MaxQuant versions that way.
The mqpar.xml that is required can be created with the MaxQuant GUI. And then the path for input files needs to be contained correctly. The path to the fasta file and the raw file. I would create a template file and then parse these filenames into the specific locations.
# mqpar.xml
...
<fastaFiles>
  <FastaFileInfo>
     <fastaFilePath>__FASTA__</fastaFilePath>
     <identifierParseRule>>([^\s]*)</identifierParseRule>
     <descriptionParseRule>>(.*)</descriptionParseRule>
     <taxonomyParseRule></taxonomyParseRule>
     <variationParseRule></variationParseRule>
     <modificationParseRule></modificationParseRule>
     <taxonomyId></taxonomyId>
  </FastaFileInfo>
</fastaFiles>
...
<filePaths>
  <string>__RAW__</string>
</filePaths>
...

The combination of maxquant=1.6.10.43 and mono=5.14 worked and a few other versions. Some versions between 1.6 and 2.x did not work at all.
Specific versions of MaxQuant are also available directly via conda and there is no need to download the zip file. Instead, you can simply install it with:
conda install -c conda-forge -c bioconda maxquant mono=5.14

Recently MaxQuant 2 has been added, and not all previous versions are available and run maxquant with:
maxquant mqpar.xml


text:
Evaluating amateur bioinformatics results
I've might have created a new algorithm for finding patterns in DNA. The technique has not been used before, and it managed to find a 33-mer with corresponding complement that exists both in E. Coli and in Salmonella. I think it's pretty neat but I don't know how to evaluate the results. In addition, I've developed this algorithm after a bioinformatics course where only the DNA sequences of E. Coli and Salmonella were provided, so I would like to test it on other sequences and have it evaluated by people who know what they are doing.
How would one get in contact with the right people for evaluating amateur bioinformatics result?

If you took a course, start by chatting with the Professor/Instructor who taught it. They will either be able to chat about things with you or put you in touch with someone locally.

It was one of those online Coursera courses. There are no real options to get in contact with the professor who created it.

What is the benefit of your approach? Is it faster than the myriad existing algorithms for pattern finding? Does it have better pattern-definition syntax? What kind of "patterns" are we talking about? What is this 33,er supposed to be? What makes it special?

I'm exploring the existance of superpermutations (of n = 4 = 33-mers) in DNA. I've created an array of all possible nucleotide superpermutations and I'm checking if sequences exist with a certain error margin. I've found a 33-mer that according to NCBI exists in 1839 unique organisms. I think this is interesting but I don't know if it means anything (with significance). The reason I'm exploring superpermutations is because they are the mother of all permutations (for which we look in DNA).

Wait, so you mean you just found some random sequence of 33nts whose only distinguishing feature is that you find it in many organisms? That doesn't seem strange, 33 is tiny, so you'll find loads of sequences that can be found all over the place at that length. Especially if you're not limiting to exact hits. Try generating some random 33nt sequences and searching for those, I bet many or even most of them will be found in hundreds of species.

Thank you. I was looking for an answer like that.

@terdon 33 is tiny, so you'll find loads of sequences I guess that depends how you look at it... with 33 nt you can make 4^33 ~ 7e19 combinations equivalent to about ~25 billion human genomes. If you pick one at random chances are it's never been seen anywhere (do I get it right?)

@dariober ah, but this isn't random at all. It's a superpermutation which I suspect will mean the sequence is relatively repetitive, making it far more likely to appear randomly in blast results. I admit I'm just speaking off the top of my head though (hence the comment instead of an answer). But yes, my assumption would be that finding a 33nt sequence present (as an imperfect match, mind you, as the OP described) in multiple genomes is not surprising per se. Finding exact matches in hundreds of genomes would be more intriguing, yes.

You need to understand the biological significance of the 33-mer locus.

I'd blast your 33-mer on NCBI and see what the biological significane is of the locus. If might be a transposon/ promoter region.
Look at the translation of the protein and if its part of a protein the location within the protein structure.
Blast the entire locus including 5' and 3' regions of the 33-mer target
Carefully assess the results for chance, e.g. conserved protein and the synonymous mutations just happened to drop out (Poison distribution).

Once you gain a greater understanding, it is likely to have been identified before and what you might contribute is a different method to obtain a known result.
In summary, you need the next level part of the course about translating a genetic correlation into a biological result/ mechanism. This isn't research grade work at present in my opinion.

Thank you for your answer, that NCBI is a really nice resource.


text:
Standalone Blast+; Automating searches, Blastn formatting, interpretation
I am struggling with using NCBI's standalone blast, particularly the blastn feature.
I have a large number of nucleotide sequences that I want to use blastn to ID. The goal is to retrieve information from the blast program and eventually import this information to excel. I want to know the capabilities of, or the easiest way, to format the output. I am running a remote search to NCBI's nucleotide database and want to filter and organize the information retrieved as much as possible. The ultimate goal is automate a way to run a remote search to NCBI's nucleotide database, and to work with the results in excel.
I am extremely new to this and struggling to find the best path. From what I can gather, many prefer using external programs, scripts, parsers, interpreters etc to organize their results. I am okay with doing this as well, but would still like to format the output of my blastn search as much as possible beforehand. It is returning a lot of data to me that I don't need.
Can I automate this in excel, and just use VBA/macros to further organize my data? It seems some people prefer to just import the output text file to a spreadsheet and work from there. I have also seen options like BioEdit, eutilities, APIs... I'm not sure which route I should invest my time into figuring out.
Thanks so much for lending your time and expertise!

You just have to change the output format using the -outfmt parameter. To have an output in a tabular format just use -outfmt 6. The resulting table have 12 columns containing the following informations queryacc subjectacc pident length mismatch gapopen querystart queryend subjectstart subjectend evalue bitscore. You may restrict the number of hits per query using -max_target_seqs 5 for example.

Hey Thomas, thanks for the reply.
I found some of these features previously but unfortunately it does not include nearly everything I want to include. The specifiers are missing some crucial information, and I would like to use more filtering options if possible. I want to somehow write something that will choose the best choice based on my criteria (probably a lot of if, else statements) to ultimately have 1 best choice for each query. I have over 30,000 queries, so I really want to minimize the amount of data that is returned to me.

I haven't found a way to have the query sequence returned in my output, nor the binomial of the subject. Since I will be having so much data returned to me, what is the best option I can use for choosing the best choice based on my criteria? Like setting thresholds and if/then statements. Also, unfortunately many of the default specifiers included for that format are not necessary for me, and many bits of info that I do want returned don't seem to have specifiers.. Its very confusing!

Hi @Sarah, I recommend you to keep the output format as a tabular separated file. You may use a -max_target_seqs 1 but it is not recommended for large searches ( https://academic.oup.com/bioinformatics/article-abstract/35/9/1613/5106166?redirectedFrom=fulltext). I think you should think about using R for example to further process your output table a write a loop to get the best hit for each and integrate if/else loops. Then you should use the  blastdbcmd to get the subject or query sequence.


text:
Outputting crunchable list of HBonds from Pymol
X-posted from the Main Stack Overflow because folks here are possibly more familiar with PyMol(?)
I have been trying to assess the strength of an interface I want to mutate using Pymol. What I am interested in is a list of contacts withing a specified interface (which I have calculated using the InterfaceResidues.py script from the PyMol wiki script library). My aim is to produce a bunch of mutants using the mutagenesis wizard and assess how the HBond network changes.
Actions>Find>Polar Contacts does the trick in visualising them, but doesn't output a file with numbers which I can then crunch to assess which mutations weaken my interface most.
I have been trying the scripts list_contacts.py and list_hb.py from Prof. Robert Campbell of Queen's University at Kingston, but with limited success (nothing gets printed or written to the file).
The scripts can be found here http://pldserver1.biochem.queensu.ca/~rlc/work/pymol/
I am going to contact Prof. Campbell directly re: scripts but, in the meantime, and in case he's too busy to deal with a Pymol noob, does anyone have any ideas of what else I could try?
I am using Pymol 2.3.1 and Python 3.7.1

Use Biopython pdb, cross check manually in Pymol a subset of the output

@Michael it's not clear to me that Biopython.PDB would be better here than PyMOL.

How do you use these script? I had only a quick look at them: running these scripts defines a new command in PyMOL. Then you need to run that command.

Covered here: https://bioinformatics.stackexchange.com/questions/12862/extract-autodocked-protein-ligand-connections-programatically/12888#12888

@marcin I need to get the list of contacts and their distances for further analysis, not just visualise them.

Thanks a lot to all who commented/replied, the author of the scripts provided a response. It was down to a syntax change in PyMol of which I wasn't aware that caused my selection to come up empty.

Can you expound upon the syntax change, ideally giving the exact line or two of code in question? That may help others running into the same problem.


text:
Is the number of contigs and scaffolds an artifact of genome assembly and rather arbitrary, or is it canonical for each species?
E.g. bowhead whale genome has X contigs and Y scaffolds, but couldn't these easily be different numbers?

Please edit your question to describe a snippet of the whale genomes in question with examples and consider the investigation more carefully. In context the answer is it is of course arbitary. How can it possibly have taxonomic value?

This is purely arbitrary and an artifact of how well the assembly turned out. Better quality assemblies with have fewer and longer contigs.

With the ideal for your example being one contig and one scaffold or each chromosome or segment of DNA.


text:
What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?
What's the difference between Barycenter Centrality and Closeness Centrality (Freeman)?
Barycenter scores are calculated as 1 / (total distance from vertex v to all other vertices) in a strongly connected network. source
Closeness Centrality (Freeman) are calculated  total distance from a node v to all the other nodes in a network. source
It looks like the same centrality .
Ps:
I want to analyze closeness for my protein–protein interaction network.

First, I'll observe that Freeman centrality as you define it is the reciprocal of the Barycenter centrality as you define. It looks like people often use the reciprocal instead of the original definition, so in that sense they do look the same.
I found the following comparison here (though it doesn't name it as Freeman):

There are 2 types of distance centrality scores, Closeness Centrality
and Barycenter Centrality. Barycenter Centrality for vertex v defined
as:
1 / (total distance from v to all other vertices)
Closeness scores are calculated using the formula 1 / (average
distance from vertex v to all other vertices) and Barycenter scores
are calculated as 1 / (total distance from vertex v to all other
vertices).

So it appears that for this package it is a difference of the reciprocal of the total vs. reciprocal of the average.
Notably, wikipedia differs here, in that it matches your source:

Freeman's closeness centrality, the total geodesic distance from a given vertex to all other vertices, is the best known example.citation

Yet another source notes the existence of a "normalized" version of Freeman's closeness centrality, which normalizes centrality of a node by the degree of the network. It also comments on whether it is the reciprocal or not:

Closeness (Freeman '79)
The graph-theoretic distance of a given node to all other nodes. The sum of the rows/columns of the geodesic distance matrix of a graph.
Simple closeness is an inverse measure of centrality: the larger the numbers, the more distant an actor is, and the less central. Should really be called "farness".
Normalized version divides the minimum "farness" possible (N-1) by "farness" to simultaneously make the range 0 to 1 and invert the measure so that larger values correspond to greater centrality (truly "closeness")

So I guess that it might depend on what tools you are using. Most of the time, it might be the same. Other times, it might be normalized in different ways with either averages or for 0-1 ranges.
Closeness centrality is probably the more standard term, so it is unlikely that anyone would criticize it. But it appears that there is sometimes some variation in what people mean by that term and whether they are using a normalized or a raw closeness centrality.
In my opinion the normalization leading to a range from 0-1 is nice because it makes the measure at least somewhat generalizable across networks. But it's up to what you need to do.


text:
MetaQuast for assembling samples from complex communities
I'm working with whole genome metagenomic samples from human skin, and I'm using MEGAHIT for assembly and MetaQuast for evaluation. However, MetaQuast requires a list of reference genomes for the analysis (if I understood correctly). When I perform the evaluation without the references, it uses BLASTN for aligning contigs to SILVA 16S rRNA database. Since I'm working with whole genome sequences, this is not an optimal choice. As expected, MetaQuast is identifying only two species in my contigs file, although I have manually checked through BLASTN that there are contigs that match other species' genomes. And this evaluation makes my assembly looks terrible. Since I'm working with samples from human skin, it wouldn't be viable to download a database of reference genomes and perform MetaQuast using a local alignment. So, what would be the best approach to perform the evaluation correctly in my case?
Looking forward to discussing it!

I would suggest first using a binning approach such as MetaBAT2 (or an equivalent tool, there are many) to separate your metagenomic assembly into bins roughly corresponding to genomes. I would then run CheckM to assess the bins' quality/completeness/purity.
If you still want to use quast, you could then run metaquast with a reference genome appropriate to each bin independently. I would suggest GTDB-tk to perform a classification and identify such an appropriate reference genome.

Thanks @MaximilianPress good way to take the question forward

Thank you @MaximilianPress, I received a similar answer at Reddit, so I believe this might be a better way to proceed indeed


text:
Download or create VCF from Human Genome Diversity Project (HGDP)
I can't seem to find HGDP data in VCF format.
Does anyone know where I can download it from?
Perhaps there a Snakemake pipeline for calling VCFs from this data?

I presume you are looking for the recently sequenced WGS data rather than the genotype data? If so, the different files can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/
Each chromsome has a seperate vcf file with a corresponding .tbi index file.
The gvcfs can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/gVCFs/
These files are more compressed and smaller in size than the regular vcfs.

Fantastic, thanks so much!

No problem. If you feel like this has answered your question, you can close it by clicking 'accept answer'.

Sure thing. When you upvote my question ;-P


text:
How do I write tests for a snakemake pipeline?
I've been playing with snakemake, which is really powerful, however, I'd like to create a 'unit test' for the whole pipeline. e.g. a way to call the pipeline that tests that all the rules work together, all the steps complete, etc.
I assumed I could find some guides about how to do this, but I don't find any clear information about how to tests pipelines.
Any advice or suggestions would be welcome.

I made a feature request here: https://github.com/snakemake/snakemake/issues/599

You can run CI on a snakemake workflow when provoding a minimal example dataset. Maybe an example from the snakemake-workflows repo will be helpful: https://github.com/snakemake-workflows/dna-seq-varlociraptor

This is what I do when I want to check whether a program runs from start to finish and produces the expected output on test data. Here I used it to test a C program but the idea can be applied also to a snakemake pipeline as in your case.
The idea is pretty simple really: use python's subprocess module to run the pipeline then check the results using the unittest module (or your favourite unit test package).
For example, save this in test.py:
#!/usr/bin/env python3

import unittest
import os
import subprocess as sp

class Tester(unittest.TestCase):

    def testPipelineMakesStuff(self):
        # Run the pipeline
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)

        # Optionally, get stdout and stderr
        stdout, stderr= p.communicate()

        # Check exits code and other expected output
        self.assertEqual(0, p.returncode)
        self.assertTrue(os.path.isfile('some.pipeline.output.txt'))

        ...Maybe read some.pipeline.output.txt and check it looks ok

    def testPipelineReturnsWithError(self):
        p= sp.Popen('snakemake [some parameters]', shell=True, stdout= sp.PIPE, stderr= sp.PIPE)
        stdout, stderr= p.communicate()
        self.assertTrue(p.returncode != 0)
        self.assertTrue('Something went wrong' in stderr.decode())

    def ...more test cases...

if __name__ == '__main__':
    unittest.main()

Run this script as:
 python test.py

and check all tests pass ok. Of course, you can embed python test.py in travis or other continuous integration system as I've done here.


text:
SNPs in LD in which populations?
I have 3 SNPs that are in LD in 'several' populations here:
https://www.ebi.ac.uk/gwas/variants/rs13078854#ld_panel_label
It's a pain to look at all (30) populations one by one to manually check if the 3 SNPs are in LD in a given population or not.
Is there a way to ask this question more systematically? e.g. in which common imputation populations are a given set of SNPs in LD (above a certain threshold)?

There is a risk of a type 2 error.
My advice would be to use Structure if it is still available given the situation. Its not great and it is old and its a GUI and I don't know whether it has been updated (I'm not exactly selling this am I?), but it will do the job of navigating through a type 2 error.
What you do is provide the number of populations a priori and it will assign a Bayesian estimate of HW equilibrium for that number of populations. You manually input 1,2,3,4,.... 30 run and wait. The point at which the joint marginal plateaus is the minimum number of populations which do not exhibit LD (thus if two of your populations exhibit LD it collapse them into one). Thus if your answer is '20' it give some confidence it is a real result. If it is like 10 (optimistically) the conclusion would be in pretty reasonable shape.
There are other ways of looking at this, but that would be my starting point.
The classic calculation is LDHat, here https://github.com/auton1/LDhat, but the last update was 2015 and in this particular circumstance I'm not sure.

I don't understand.

Sorry if this has not been helpful. I would therefore suggest a data pipeline around what appears the current standard which is plink https://www.cog-genomics.org/plink/1.9/ld . To understand LD stats a background in population genetics and stats is really helpful. The answer above is as good as I could do without writing a data pipeline.

We can chat about it, but the solution I proposed was trying to over the type 2 error in the data and I'd need to explain that in depth because its fundamental to how LD works. I'm not up to date with plink, I left LD prior to plink, but I do understand how the numbers work in LD.

How would you like to chat?

I forget how you do invites and stuff, but I think you have to visit the bioinformatics chat and that activates an invite.


text:
Is it valid VCF not to 'squash' positions with more than one ALT allele?
I'm seeing output from PBWT that looks like this:
3   112841  rs78923776  C   G   .   PASS    RefPanelAF=0.000153988;AN=2;AC=0;INFO=1 GT:ADS:DS:GP    0|0:0,0:0:1,0,0
3   112841  rs78923776  C   T   .   PASS    RefPanelAF=0.078534;AN=2;AC=1;INFO=1    GT:ADS:DS:GP    0|1:0,1:1:0,1,0

The first line says that my genome is CC at this position, but the second line (for the second alt allele) says that my genome is CT at this position. OK, the first line couldn't say this, so it calls me as REF/REF, but this call has to be interpreted in the context of the second call, and can't be taken at face value.
That's why I wonder if this is valid VCF or if it's actually breaking the specification.
The line should be 'squashed' down to the following (I think):
3   112841  rs78923776  C   G,T .   PASS    RefPanelAF=0.000153988,0.078534;AN=3;AC=0,1;INFO=1  GT:ADS:DS:GP    0|2:0,0:0:1,0,0

Many pieces of software can only work with biallelic diploid calls, so they have to split up multiallelic genotypes into separate lines.

Thanks. Is there a standard way to work out the genotype call?

It is a valid multi-allelic entry, and most variant callers produce such entries. However, comparison between two such gets more complicated. With biallelic entries, 2 CHR-POS-REF-ALT entries either match or they don't. With multi-allelic entries, they can match completely, partially or not at all, which means you also need to be able to specify how much of a match is a partial match etc. This is a hindrance in annotation processes that annotate using the CHR-POS-REF-ALT fields.
This also introduces cases where an indel can become a multi-allelic entry with a nearby SNV like so:
chr1    100    GA    G
chr1    101    A    AA

can now be represented as
chr1    100    GA    G,AA

which can get confusing. Operations are always more predictable when VCF has been decomposed and left-aligned.

Good to know this is valid. I'm having trouble interpreting some of the genotypes I'm seeing. Is there a standard way to do this?

Can you give me an example of the difficulty you're facing?

Hi @Ram, there are several genotype call combinations that I don't know how to 'decode'. I've written up a report here:
https://github.com/richarddurbin/pbwt/issues/52#issuecomment-1064125959
The simplest way I can describe it is to say that, if the call is split across two lines (for a tri-allelic variant) the REF part of the call should be the same on both lines, no?

I glanced through your thread there and I think you're in the gray area between biology and bioinformatics. VCF entries assume reference sequence is reality, so yes, a multi-allelic site, when split, will have the same REF bases. Assuming that the mutated site is still diploid, a tri-allelic variant will have GT 1/2 where 1 is ALT-1 and 2-ALT-2 (so no REF base)

Right, I'd put that in VCF as something like ./1 and ./1 (assuming I had to split alternative alleles out into separate bi-allelic lines for some reason).

Yes, that's what decomposers do - they give you 1/. GT entries.


text:
Given a 'random' rsID, how can I find the closest (LD) SNP in a given set of SNPs?
I'm assuming that this is such a common question that I don't have to do my own imputation...
I have a SNP from a paper that is associated with a specific trait... Actually I have 8 of them:
https://www.ebi.ac.uk/gwas/publications/33307546
However, when I search my SNP data from 23andMe (v3 NOT v5), only one of them is in the set of SNPs used by that chip, namely rs2109069.
I'm wondering how I find the closest SNP in my v3 data for the remaining 7 SNPs missing in that data.
Is there a quick way to find a 'tag' SNP in v3 for a given 'random' SNP?
I'm aware that a tag SNP will only be as good (on average) as the LD between it and the given SNP.

What is a "closest SNP"? You mean in distance? Why would that be relevant? Don't you want to look for pathogenic or at least LP or VUS variants? Why would the nearby SNP be of interest?

Sorry, thought it was clear. SNP x is associated with a trait but isn't on the SNP-Chip. SNP y is in LD with X and is on the chip.

So you want to try and guess at the variants your sample carries based on their linkage disequilibrium values? And is this using the raw, unfiltered data of 23&me? (if the answer is yes, please don't do this, their raw data are not supposed to be used that way and are not reliable at all).

I don't understand... you mean I cant run 23andMe v3 data through imputation?

No, I mean you shouldn't trust 23&me data beyond the filtered results they show you. If you have downloaded the full, unfiltered dataset, that will contain loads of false positive calls that should not be considered. Imputation won't be a problem, I just find it a bit odd that you would use what is essentially a statistical approach (using LD values as a proxy to guess what variants are really present) when looking at a single, specific sample. Is getting proper sequencing data so you can actually see what you have in the loci of interest not an option?

Yes it's not an option. Please can you provide a reference for the fact that 23andMe 'unfiltered' data contains a lot of false positives?
I had understood that SNP-Chip accuracy was very high when validated by WGS techniques.

I don't have a reference, but I have looked at the unfiltered results of four different (healthy) samples and all of them had variants that, if true, would have made these healthy individuals very, very sick or, more likely, dead. Variant calls need to be filtered and the 23&me unfiltered results are just a list of variants with no information on call quality, coverage, fastq quality or anything that can be used to sort false positives from true calls. That's why they only show the filtered results.

You're talking about sequencing based variation calls. I'm talking about SNP-Chip based calls. As far as I know, 23andMe isn't doing any sequencing commercially.

All I can tell you is that in all 4 cases I have checked, the unfiltered data, the dataset you can download but is not shown in their web interface, was riddled with obviously unreliable calls.

That makes no sense. For example, I ran my own 23andMe data through GenomePrep (results here: https://supfam.mrc-lmb.cam.ac.uk/GenomePrep/results.html?id=48318966999) and you can see that 80% of the calls match the reference. Sure 20% don't, but I wouldn't be me if they didn't.

As you wish. I can only share my experience based on what I have seen and my ~7 years of working in NGS variant calling in a clinical setting. I would never trust unfiltered 23&me results. Their filtered stuff is fine, but unfiltered variant calls are generally not reliable.

That's true when it comes to GBS data, but as I said, 23andMe DON'T do GBS, they use SNP-Chips, which are generally regarded as very reliable.

And I'm sure their filtered data are indeed reliable. Not so their unfiltered data, at least not in my experience. There's a reason they don't show the unfiltered results, after all. Here's a paper that seems to confirm my experience (I only glanced at the abstract): https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7879796/: "SNP chips are extremely unreliable for genotyping very rare pathogenic variants and should not be used to guide health decisions without validation."

This is weird... You don't seem to be understanding what I'm saying, but you clearly are very good in the field, so I'm confused. Perhaps we should have a call, because we seem to be talking over each other. e.g. of course a SNP chip won't see a rare variant... it's totally beside the point... Somehow I think we're talking about different things.

Heh, I'm getting the same impression :)  Feel free to ping me in the site's chat room. Might be easier there.

Also, while I do have experience with NGS in general and NGS in a human clinical setting in particular, I do not have any to speak of with SNP chips. I just know that the unfiltered 23&me data I have seen have consistently been crap. I may very well be wrong or may have misunderstood what you're doing entirely.


text:
Is SNPRelate the goto package for PCA of SNP data in R?
I'm trying to evaluate a method that I'm not familiar with. They used SNPRelate. I'm wondering if this is the best / only choice... what other tools would you recommend?
Many thanks,

Does it have to be in R? If not, the standard things to make a PCA from genotype matrix is plink or eigenstrat.

That said, unless you are doing fancy stuff like PCA projection - PCA is just PCA and whilst methods might differ in terms of computational compexity, the resulting eigenvectors will be mostly very similar.

SNPRelate claims nearly 10x faster than plink... but I don't know.


text:
Convert SNP data from GEO into PED and MAP files?
I can easily grab SNP-Chip data from NCBI's GEO database for a given Illumina SNP-Chip array, e.g. from here:
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52147
GEO provides data as SOFT formatted family file(s), MINiML formatted family file(s) or Series Matrix File(s).
I'm looking for a way to convert these into PED and MAP files for use in PLINK...
Or does PLINK read these formats natively?

I see. Maybe PLINK cannot handle this. I see on GEO that series downloads are suitable for opening in spreadsheet.

I don't want to manually play with 3000 sheets...

Sadly the 'contact us' email at GEO didn't help at all.
I ended up probably re-inventing the wheel with this script:
https://gist.github.com/CholoTook/60968e3ab6d90cb8fd19be55a25592f1
YMMV


text:
Inflated p-values in quantitative trait analysis
I am performing a quantitative trait association between the expression of one gene and ~400,000 methylation values. First, both variables are rank inverse normal transformed, adjusted for confounders (technical and for family structure) in a mixed-effects model and we use the residuals for the association. Then we do a linear regression between the residuals of expression and each methylation site, adding other four confounders in the model (age, BMI, etc.). The sample size is of ~500. For the association we are using MatrixEQTL.
When we check the p-values, they are clearly inflated. Is it something that I have to be concerned about? Did I probably miss something in the analysis? Thanks in advance.

it's a long chunk of data processing... so normally it looks like this when you are using the wrong distribution... what is the test used by MatrixEQTL ? One way to check whether that is performing ok is to permute your residuals and run MatrixEQTL again.. If it doesn't follow the diagonal, it means your data is not suitable for this.. amazing package

It performs a linear regression (least squares model) and a t-test for the significance of the coefficients.

try permuting your phenotype and see whether this qqplot still holds?

Thanks @StupidWolf

@MaximilianPress explained it better than me lol... so if you permute does this inflation go away? Honestly, I have always done it from scratch.. this chunk you mentioned above, i supposed involves some god zillion packages? Hard to tell where it goes wrong

It cannot be that your whole genome is associated with that 1 gene expression??

Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.
If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.
Some random thoughts:

What does the distribution of your methylation data look like after normalization?
Your observed distribution is basically never on the expected line, which argues strongly for confounding.
How are you adjusting for non-independence of methylation types? Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model, or at least fit 10-20 principal components to account for population structure.
How are you representing your methylation status? If it's binary you might look at this paper.

Thanks @MaximilianPress. Actually I'm don't think we have other confounding variables, we are correcting for all the possible ones, including family structure, and all data come from the same population.
After the adjustment for residuals of methylation, data is normally distributed. It is not binary, but the probability of being methylated. I'll check the papers you mention.

@SergioVillicana  You honestly can never know whether you have corrected for all confounding variables. Family structure is ok but underlying genetic relatedness, e.g. cryptic relatedness, is much more difficult to detect. Additionally, I realized that the Arabidopsis paper (mentioning the epigenetic similarity matrix) might not be appropriate, I got confused thinking about methylation as response vs. explanatory. I'm editing the answer to remove that paper.


text:
Odds ratio and enrichment of SNPs in gene regions?
I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).
A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:

Take a single random sample of all SNPs with size 10n (n is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).
Do a bootstrapping, taking multiple random samples (excluding qtl-SNPs), each of size n, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a p-value.

Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?
Thanks.

also posted in https://stats.stackexchange.com/questions/485997/odds-ratio-and-enrichment-of-snps-in-gene-regions. Avoid posting the same question twice. I would suggest thinking about the nature of  your question which forum is more suitable

@StupidWolf the MAF distribution of the full 7M and the qtl-SNPs is different. That's why I don't want to consider the full panel and instead use a more comparable 'background' of SNPs to have a meaningful estimation of the enrichment.

The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.
My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.
This is what I think might work

Calculate your true OR, qtl snps versus all others for say TFBS
Divide your snps according to MAF into 10 bins / stratas
Within each bin, we permute the label whether the snp is a qtl
Calculate OR of selection / others for TFBS
Repeat 1-3 N times to collect the permutated OR
p value would be number of (permutated OR > observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

We have the MAF for all snps, with those in qtl being different from others and put them in bins.
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

We sort the snps according to MAF to make it easier:
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

And we create the bin:
bin_MAF = cut(snp_MAF,10)

Lastly a function to calculate OR:
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

So our observed OR (in log) is:
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

Now we collect the OR under permutation:
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)

Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation


text:
extracted counts matrix from Tximeta generalized experiment with no linked transcriptome are normalized?
I used the Tximeta vignette to construct a summarised experiment from quant files (salmon, index prepared with NCBI transcriptome).
1)Tximeta return me a non-ranged se (which I am not sure about the real implications)
because of this...
2) I used makeLinkedTxome to link the transcriptome and the gff (local files used in the index creation).
3)run tximeta on the colData to create the summarised experiment se
4) seqinfo(se) returns me a NA populated table
5) I summarise to genes (gse)
6) I extracted the gene counts counts <- assay(gse)
My questions here are;
a) Are these counts normalised?
b) If not, how can I normalise them? TMM for example? (CalcNormFactors will not work here right?)
Is there any better way to obtain a normalised count matrix from tximeta? (maybe I don't need to link the transcriptome...)
Any help would be appreciated!
Thank you! and regards,

Thank you! so I can use these `counts'  with the pipeline. Perfect! thanks!!!


text:
apply TMM on counts imported from salmon using Tximeta
I used Tximeta to import a summarisedExperiment from the salmon output (used with genocide transcriptome v34). I need to produce 4 matrix of counts:
- tx counts in TPM
- gene counts in TPM
- tx counts normalised
- gene counts normalised
The first two I used assay(se, "counts") and assay(gse, "counts") respectively.
Mu problem is to apply the normalisation to them. I used deseq2 to estimate the size factors with a fake design because there are no groups yet in the dataset so I used ~gender to been able to generate a adds object to apply the EstimateSizeFactor, but when I try to see them returns me 'NULL'
> dds_genes <-DESeqDataSet(gse, design = ~ gender)
> dss_genes <- DESeq2::estimateSizeFactors(dss_genes)

> DESeq2::sizeFactors(dss_genes)
[1] NULL

I get count values like this after using these commands (which I am not really sure they are really normalised;
normalized_counts <- DESeq2::counts(dss_genes, normalized=T)
normalized_counts[1:10, 1:10]

ENSG00000000003.15  654.92903  771.34051  500.1863  560.30530 1144.8295  938.40446 1173.208233  896.28964  388.1482   501.20545
ENSG00000000005.6     0.00000    0.00000    0.0000    0.00000    0.0000    0.00000    1.344747    0.00000    0.0000     0.00000
ENSG00000000419.12  416.59305  399.05434  675.0911  285.64627  551.8253  424.59423  461.448691  679.08099  365.8343   461.66276
ENSG00000000457.14  242.52226  244.91792  279.9109  414.24494  262.8790  354.50834  451.671501  427.08499  407.7070   329.50667

I don't know how to proceed to just get a TMM normalised values from tximeta object. That would be my preferred option. I guess deseq2 can also provide me with normalised data.
Maybe there are other ways?
Thank you for the help!!
Regards


text:
method to compare the similarity between proteomic profiles
I have a problem finding a methodology to compare proteomic profiles against another. The experiment is the following:

Experiment 1. Ten pdx tumour models with their proteomic profiles were done in one Research Institute with DIA-MS. All of them are the same tumour from different patients. Differential analysis was done with Limma.

Experiment 2. Same tumour, but this time entirely on humans (not pdx) and also controls were included (SWATH-MS in another research institute). Differential analysis was done by t-test.

The biological question is which of the 10 pdx models resembles better with the human one.
I tried different ways of normalizing together the raw data but there is a huge difference in the levels of abundance between both experiments. Also, the pdx ones are really similar between them.
In the end, I did a one vs rest contrast on the pdx conditions to get the unique markers for each model and use them as genesets (positive and negative separated) to try to see if the upregulated and downregulated genes (tumour vs control on experiment 2) are some kind of enriched (pre-ranked GSEA with rank as -log10 p-value*log2FC).
Do you think this can be a solution? Is there a more simple way to do this assessment?
Thank you!!


text:
Receptor-Ligand R-package for visualization
I was wondering if there was any R-package for receptor-ligand interaction similar to SingleCellSignal (https://github.com/SCA-IRCM/SingleCellSignalR_v1/blob/master/SingleCellSignalR/vignettes/UsersGuide.Rmd)
+++ if there are vignettes and/or the graphs can be customized.
+++ if it can be linked to Seurat
Note: I have tried iTALK, but I keep getting an error that seems to be a common problem and have left a query.
Thanks

Two nice tools for receptor-ligand interaction analysis are:

CellPhoneDB by the Teichmann lab, although this is not an R package...
NicheNet by the Saeys group

Note that these tools are mainly meant for analysis of receptor-ligand interactions, although I think NicheNet also has some visualization options. CellPhoneDB needs to be run from a shell, but the results can easily be imported into R for visualization. You might check the papers on these 2 methods for some ideas of visualization. Combining Seurat's visualization options and ggplot2 should get you already pretty far.


text:
plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
I have used the empirical RUVg method (from RUVseq) to estimate the unwanted variation of my dataset (consisting of several public datasets analysed together, with controls and case samples in different datasets.
Afterwards, I used EdgeR to perform my DE analysis with 2 RUV covariates:
designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

I got my results after making all the contrast that I needed but I would like to plot the expression of some genes of interest. One plot per gene with all the conditions in condition4 has 9 categories (expression level or norm counts on the Y axis).
My doubt is which counts to use. I have the normalised ones from RUVg, but it states they should only use them for exploratory analysis (maybe this is the case??). These counts as I understand are the residuals of the model, so not really sure what to do.
Another option is to use the normalised counts produced by EdgeR, but these ones doesn't account for the unwanted variation (W_1, W_2 covariates) (right???)
Any advice or nice solution?? :=)
Thanks in advance!!

You could correct your logCPMs (edgeR::cpm(y, log=TRUE)) with the covariates from RUVseq using limma and then use these for plotting:
set2 <- RUVg(...)

designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

# correction with limma
logcpm <- cpm(y, log=TRUE)
design <- model.matrix(~0+condition4, pData(set2))
corrected <- removeBatchEffect(logcpm, design=design, covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])

Thank you @alderaan_shot_first. Can I ask why you used a "-" in the model.matrix? Also, after reading your answer I have been checking the removeBatchEffect documentation and when you put the covariates argument you are setting all the phenoData of set2... Don't would be covariates=as.matrix(pData(set2)[,c(W_1,"W_2"]) ?. Thank you!!!

was a typo, sorry. As for the matrix, yes it is fine with the code you suggest. feel free to make an edit.


text:
about the scaling after normalization and transformation in RNA-seq data
When we use count data in RNA-seq analysis, we usually use normalization and sometimes vst, rlog transformation (DESeq2)or log2(CMP+4) transformation (edgeR) to perform K-means clustering. Can we use scaling after normalization and transformation?
Or Should we use scaling only when we just use normalization without transformation?


text:
about the normalization of RNAseq data for calculating distance for unsupervised learning
I have been working on clustering using RNAseq data. To compute distance, what kinds of normalization is optimal?
Can we use normalization using relative log normalization? I understand this should be used for differential gene expression analysis, not for unsupervised learning in which we do not use information of disease condition.
If you can show good papers, please share.

Hi it would be useful to describe the disease in question and some idea of the data. "Normalisation" is often used in parametric statistics to convert data to a normal distribution as a pre-requisite of the analysis, or linear modelling. Transformation is a "hit and miss" business in general. Unsupervised learning ... it depends what you mean, but you wouldn't do this for PCA or tSNE. More details please.

Typically some sort of variance stabilizing transform is used before clustering. Popular options are regularized log transformation or a vst transform, which are available in DESeq2. Note that these are NOT used for performing differential expression, just things like clustering and PCA. For differential expression one would use TMM, RLE or similar.

Thank you for your comment. Yes, I would like to calculate the distance to perform clustering, not differential expression.

So if we want to clustering after calculating distance, you mean we should use regularized log transformation or a vst transform. And when we want to perform differential gene expression analysis, we should use TMM, RLE. So we need to change normalization methods in one paper depending on analytic plan, right??

Yes, though regularized log and such are transforms and not normalizations.

So for the analytic pipeline to compute the distance, we should not use both normalization and transformation ,right?  raw count (without transformation)→vst or other transformation.
https://towardsdatascience.com/pitfalls-of-data-normalization-bf05d65f1f4c

Please see and follow the DEseq2 tutorial, this is covered there.


text:
about two group comparison of three group data in DESeq2 package
When we have three groups samples (A,B,C) with 25000 genes and the main interest is A vs B, should we limit samples only A and B for normalization to perform DEG analysis? Or better to include all samples (A,B, and C) for normalization and specify the comparison group within three groups?
res <- results(dds, contrast=c("tissu_type","A","B"))

I guess it makes more sense to include A and B. Imagine you are writing a paper, and it simply does not make sense if you include C, which you never talk about in the paper. But there is no absolutely correct answer for this. I also believe that the results are going to be very similar. If the results using two methods are so different, then there might be something wrong.

Thank you so much for your all comments!

As explained on the comments, having some other samples might distort your estimations, add an outlier or something alike.
The advantage of having all the three (or all) groups of samples before testing differences, because usually methods like limma, DESeq2 and others uses the information of variation of the genes across all samples to update their prior probabilities and the variation of the genes on the statistical model used to compare between groups.
If the method and your model are robust use all the three groups (generally the methods used are robust, but sometimes the models are not). By robust I mean that it doesn't miss any know (and available) source of variation. Otherwise the comparison might be comparing only the samples not the effect of the condition on the samples. If you only compare the samples it would be equivalent to use a t.test to the limma methods (not completely sure on the DESeq2 case). Or explained differently, you wouldn't take advantage of the statistical improvements implemented on these methods
Usually the reasons why these groups are sequenced on the same batch is because they are related (or to add some known ground truth), so you can release it with the data when communicating your results but you may omit that you used this third group on the methods section.

Plot out the PCA of your samples.  If, say, PC1 splits the C's away from the others, and PC1 is more than, say, 50% of the variance, or if the C samples are a lot more spread out than the A and B samples, keep it apart.  Otherwise, do them all together.
If the C's aren't too different from A and B, the math will benefit from including them.


text:
Available Protein sequence alignment dataset and HMM model
It may better to move the question here.
I am new to biology and I find my algorithm may be used in the Protein sequence alignment, since it is a henced HMM model. I find that people use HMM to generate noisy copies of the consensus sequence of different lengths. There is a figure show the process:

I am familiar with the machine learning or graphical model. To build the HMM, we need to set the model parameter, which we learn from the training data in ML. My question is whether the train&testing works for this question? If so where can I find data to training this HMM model?
It seems Professor Richard Durbin may release some dataset. But how can I find available datasets? I lost in the bio-vocabulary and fail to find datasets. Also I am wondering whether or not this topic is very minor in this community?
Update:
Since I may abuse the vocabulary, it may be better to post the original word I read.

Why don't you just build an MSA (with MAFFT/MUSCLE/CLUSTAL-OMEGA) and then build a pHMM with HMMER3? Or if you want already existing HMMs, checkout PFAM


text:
calculateQCMetrics defunct, how to calculate the quality metrics by perCellQCMetrics
I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use calculateQCMetric function to calculate the quality metrics, but I am getting an error message that calculateQCMetric is a deprecated (isSpike function is also deprecated). I want to calculate the quality metrics. Can anyone provide code for getting QCMetric step? I attached the previous defunct code below:
umi <- calculateQCMetrics(
   umi,
   feature_controls = list(
       ERCC = isSpike(umi, "ERCC"),
        MT = isSpike(umi, "MT")
   )

And I use StupidWolf's code from this question isSpike function in SingleCellExperiment package is deprecated? to instead isSpike function. I also attached below:
Spikein_names = grep("^ERCC-",rownames(molecules),value=TRUE)
SpikeIn = molecules[Spikein_names,]

g = genes(EnsDb.Hsapiens.v86)
MTgene_names = g[seqnames(g)=="MT"]$gene_id
MTgenes = molecules[rownames(molecules) %in% MTgene_names,]
Have a vector of genes to keep:

keep = setdiff(rownames(molecules),c(SpikeIn,MTgenes))
Then make the single cell object:

sce <- SingleCellExperiment(
assays=list(counts=as.matrix(molecules[keep,])),colData=annotation)
Use altExp to slot in the stuff:

altExp(sce, "spike-in") <- SummarizedExperiment(SpikeIn)
altExp(sce, "MTgenes") <- SummarizedExperiment(MTgenes)

I try the below code, which did not work.
> perCellQCMetrics(sce,
 subsets = list(Spikein_names,MTgene_names),
 flatten = TRUE,
 exprs_values = "counts"，
 use_altexps = TRUE)

Thanks in advance.

I'm also following the Sanger tutorial. I've tried the code suggested by mattcwh but get the following error: Error in assay(altExp(x, use_altexps[i]), exprs_values) : 'assay(<SummarizedExperiment>, i="character", ...)' invalid subscript 'i' 'counts' not in names(assays(<SummarizedExperiment>)) It might have something to do with logcounts, at least based on this post (?) : https://github.com/LTLA/scater/issues/3

Ok I have figured it out and I am sharing this in case anyone can get help. I used
sce <-addPerCellQC(sce)

which calculated the QC for all the altexps also and adds them to the colData.

It would be helpful if you can provide the error message.
Guessing from your codes, it's possible that the alternative experiments weren't set correctly. My suggestion:
altExp(sce, "spike-in") <- SummarizedExperiment(assays=list(counts=as.matrix(SpikeIn)))
altExp(sce, "MTgenes") <- SummarizedExperiment(assays=list(counts=as.matrix(MTgenes)))

If you're using altExp, the subsets argument in perCellQCMetrics is unnecessary in this case.
perCellQCMetrics(sce,
exprs_values="counts",
use_altexps=TRUE)


text:
Error in snakemake when running conda envs in container
I am facing an issue when running snakemake with a container and conda.
Let me explain better.
I am developing a workflow (locally) in snakemake, where some rules use a conda env (defined in /envs/env.yaml) and other rules use wrappers (so they also are based on conda envs).
An example here:
rule clean_spike:
    input:
        sample_ref="{}results/bam/{{id}}.tmp.bam".format(outdir),
        sample_spike="{}results/bam_spike/{{id}}_spike.tmp.bam".format(outdir),
        sample_ref_index="{}results/bam/{{id}}.tmp.bam.bai".format(outdir),
        sample_spike_index="{}results/bam_spike/{{id}}_spike.tmp.bam.bai".format(outdir),
    output:
        sample_ref="{}results/bam/{{id}}.clean.bam".format(outdir),
        sample_spike="{}results/bam_spike/{{id}}_spike.clean.bam".format(outdir),
    conda:
        "../envs/pysam.yaml"
    log:
        "{}results/logs/spike/{{id}}.removeSpikeDups".format(outdir),
    script:
        "../scripts/remove_spikeDups.py"

Locally, the pipeline runs smoothly.
Now, I wanted to test the pipeline on another machine (server from the institute), and so I run the command
snakemake --containerize > Dockerfile
I got a dockerfile and  I generated an image that was uploaded on dockerhub here.
At this point, I copied the workflow on the other machine, I added to the Snakefile the line:
containerized: "docker://davidebrex/chiprxsnakemake:latest"
and tried to run the pipeline with: snakemake -c 10 --use-conda --use-singularity.
Ideally, snakemake should run the conda envs within the docker container for each rule
The problem is that i get an error related to the fact that conda is missing..
I checked the docker image and I can run conda and see the installed envs.
Any suggestions about what is going on?
Thank you!
This is the error I get:
/usr/bin/bash: conda: command not found
Traceback (most recent call last):
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/__init__.py", line 792, in snakemake
    success = workflow.execute(
              ^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 1246, in execute
    raise e
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 1242, in execute
    success = self.scheduler.schedule()
              ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/scheduler.py", line 636, in schedule
    self.run(runjobs)
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/scheduler.py", line 685, in run
    executor.run_jobs(
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 158, in run_jobs
    self.run(
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 582, in run
    future = self.run_single_job(job)
             ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 635, in run_single_job
    self.cached_or_run, job, run_wrapper, *self.job_args_and_prepare(job)
                                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/executors/__init__.py", line 623, in job_args_and_prepare
    self.workflow.conda_base_path,
    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/workflow.py", line 456, in conda_base_path
    return Conda().prefix_path
           ^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/deployment/conda.py", line 653, in __init__
    shell.check_output(self._get_cmd("conda info --json"), text=True)
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/site-packages/snakemake/shell.py", line 61, in check_output
    return sp.check_output(cmd, shell=True, executable=executable, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/subprocess.py", line 466, in check_output
    return run(*popenargs, stdout=PIPE, timeout=timeout, check=True,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/davide.bressan-1/mambaforge/envs/new_snakemake/lib/python3.12/subprocess.py", line 571, in run
    raise CalledProcessError(retcode, process.args,
subprocess.CalledProcessError: Command 'conda info --json' returned non-zero exit status 127.

I tried to change the snakemake command and remove --use-conda but then I get an error related to the fact that the tools are missing.
I checked out other questions on stackoverflow but they were not helpful in my case (like this one: here).

I have zero experience with --containerize, but "I checked the docker image and I can run conda" versus "conda: command not found" when running the workflow really smells like a conda shell setup issue to to me, like people run into when trying to call conda in a script versus interactively.  (No idea where to point you from there, though!)

Thank you. At the end I changed snakemake version and it worked..quite weird

Actually now i still get the same error, I also opened an issue on snakemake github, I am clueless


text:
Plink 1.9 equivalent of --score sum in Plink 2.0
I have previously performing Allelic scoring in Plink 1.9 using:
plink \
--score sum sumstats.txt\

However, it seems this has changed in Plink 2.0. What is the equivalent function for Allelic scoring?

plink2 --score cols=+scoresums

See https://www.cog-genomics.org/plink/2.0/general_usage#colset and https://www.cog-genomics.org/plink/2.0/score for details.

Hi a bit of explanation would be helpful

I don't have enough reputation for a comment but the place where the option to add scoresums is listed is here https://www.cog-genomics.org/plink/2.0/formats#sscore


text:
Converting aligned fasta to plink ped/bed
I have an alignment of multiple sequences in a FASTA file (output from MAFFT), for which I would like to simulate a phenotype using plink, but for that I need to have my alignment in a PED file or plink BED. There's plenty online about converting in the other direction, but the only mention I found of converting FASTA to BED is about aligning a raw FASTQ file for use in plink - my data is already aligned! How can I do this conversion?
Note that How to convert FASTA to BED is about the UCSC BED format, not the plink binary PED.

FASTA file contains sequences, whereas PED file contains variants. Usually, you need to do variant-calling on sequencing file to get variant file. I do not know that is contained in your FASTA file.

I have the same challenge, I have aligned sequences but not vcfs, bams or fastqs. Any ideas?

An alignment can be the result of two slightly different analyses:
There is multiple sequence alignment (which is what you get from MAFFT) where sequences are aligned so that similar regions are on top of each other. This may require introducing indels (insertions / deletions) if a particular region is absent in some of the sequences.
Then there is alignment to a reference sequence. This is usually how you get variant data and most of the tools deal with this use case. Therefore, in the link you provided the solution includes aligning the data to a reference sequence. Usually the refernce sequence is much larger (e.g. the human genome) and after alignment the result is a file (.bam) that tells you where the query sequences match in the reference.
I did a quick search for ways of converting a multiple sequence alignment (MSA) to VCF. There is a tool called msa2vcf in the Jvarkit collection of utilities that can do this. The example is for the CLUSTAW format but FASTA is accepted as well.
You shoud check that indels were treated correctly because these are most likely to cause trouble.
Then you can simply convert the VCF to PLINK format  using PLINK.


text:
MEME's motifs have an E-vaue too high
I'm currently learning how to use MEME-Suite (web tool) to find a specific known transcriptor binding site in the promoter region of 5 putative bacterial genes of interest, but the results obtained with just the five promtoer sequences have an E-value too high. What can I do to improve these results? Can I use as input some other reported sequences that have the motif I'm looking for without introducing any bias?
I used the option ZOOPS, with a background model of first order, and I limited the search to palindromes only, since I'm looking specifically for a 20 bp sequence with dyad symmetry. The input sequences are the promoter regions of each gene (200 pb upstream of the start codon).


text:
GWAS Rooted PCA analysis problem
I'm fairly new to plink software and wanted to get some additional practice after doing several tutorials. I obtained the data from this paper (I'm not using this paper's methods) to do some QC with subsequent analysis because they mostly use SNPstats.
Some background on what I did leading up to PCA:
My QC steps involved:
Filtering for low MAF
Filtering for high missing values on the individual and SNP level
Excluding outliers in heterozygosity (based on F statistic)
Excluding controls that deviated from Hardy-Weinberg
Finding related individuals and excluding members from the pairs that had the lower call rate
After these steps, I had 657649 variants remaining, and no individuals had been filtered out from the original total (1401)
I then used a 1000 genomes dataset that I had already downloaded (unsure of which phase, but it contained 690 individuals) and merged it with my data (did similar QC steps on the 1000 genomes data, and also used the --flip flag to fix strand issues).
I did a rooted PCA analysis (with a LD threshold of 0.2, which included 112795 variants post-filter, but I also tried PCA without the filter with similar results) to look at population substructure and got this plot (in R):
Where the pink "OWN' labels in the plot are my data points. Does anyone have any idea why my data wouldn't be lining up with any of the ethnic groups from 1000 genomes? The data is said by the authors to be "pre-filtered," which I'm assuming means fairly homogenous, so I wasn't expecting them to be in different ethnic groups, but they should be lining up with one of them, right?

It would be helpful if you attach your resulting plot here.

Agreed 100% with @PhoemixMu


text:
nanopore QC measures on fastq file
I am working with virus nanopore sequencing data. I have to remove host genome from the reads. I want to run the QC on the fastq file after host genome subtraction. I have used nanoplot. It is slow. I dont have summary.txt. Is there any nanopore read specific qc package which can directly take nanopore fastq files.

Oh yeh, NanoPlot can be pretty slow when you only have fastq data. I'm happy to look at contributions on how to make it faster. Do you have a machine with multiple cores? Because if your data is split across multiple fastq files then NanoPlot can read those in parallel.

I have 6 core with 12 threads. How can I split the fastq files and run them in parallel?

Splitting can be done e.g. with https://fastqsplitter.readthedocs.io/en/stable/, or the Unix split tool.


text:
Regarding RNA seq data analysis and building coexpression network
I have some questions regarding RNA seq analysis if you can suggest anything it will help me a lot.

I am currently normalizing RNA seq data for comparing genes expression within and between samples. Which normalization method would you recommend for this type of analysis? FPKM, TPM, TMM? Also I want to make a heatmap to see genes expressed in different conditions. Do you think transforming normalized data (like log2, z-score) is a good idea for this?
Also I want to build a co expression network, I am just wondering if normalization like FPKM, TPM, TMM  has any influence on building a coexpression network?
Another thing I want to do is use Pearson correlation but I am confused if it will only find linear relationships among the normally distributed data.  But normalization methods do not assume that counts to be normally distributed. So, Is it a bad idea to find pairwise coexpression of genes using pearson correlation? If so, which method would you recommend is reliable for building coexpression networks?

Please help me with this, I will be very grateful to you
Thanks in advance


text:
median normalization for proteomics
I am using the data from a proteomics study were the data was log2 transformed and then a median normalization was applied. The data was normalized by groups of conditions (normal, mutant), not for all the arrays at the same time, so the median for the normal group is slightly higher than for the mutant group.
I am wondering if this is correct or should I normalized the data en general, I mean taking into account all the arrays and normalizing all at the same time so the 2 groups have the same median.
Regards

The kind of normalization depends on what you what to explore. So, there is no absolute answer to this. Different normalizations highlight different aspects of your dataset. There is a nice paper for a higher throughput quantitative study that I would recommend reading:
Nusinow, David P., John Szpyt, Mahmoud Ghandi, Christopher M. Rose, E. Robert McDonald 3rd, Marian Kalocsay, Judit Jané-Valbuena, et al. 2020. “Quantitative Proteomics of the Cancer Cell Line Encyclopedia.” Cell 180 (2): 387–402.e16.
What you want to aim for is a normalization that removes batch effects and that clusters together your quality control samples that you know are very similar to each other.


text:
Non-exome variants called from whole-exome sequencing data
I'm working on rather standard whole-exome sequencing data and treat it the same as whole-genome sequencing data (aligning it to the full GRCh38 reference assembly and calling variants with no exome-related parameters).
I recently noticed that my variant calling results contain numerous variants mapped to a non-exomic region, incidentally, some of them were significantly associated with the phenotype I study.
My questions are as follows:

Are the non-exome variants in my exome data invalid? or could it be that the exome kit captured more than it was designed to capture?
Should I treat raw exome sequencing data any different than genome sequencing data, and if so, what parts of my pipeline should be modified to accommodate exome data?


text:
Analyzing transcripts for non-coding RNA
From a RNA seq experiment we look for the counts of transcripts of tRNA but we are also interested of other types of RNA we can find. Because of the protocol used for the library prep we can't analyse miRNA but we can analyse ncRNA.
I never analyzed transcripts for ncRNA and I have some questions:

The name of the transcripts are like ENST00000611544.1, ENST00000387069.1, since they are non coding, are they linked to a gene? Should I rename the name of each ncRNA with a gene name or working with these long names is ok?
When analyzing transcripts of RNA I did enrichment analysis to see what pathways are more active. But when analyzing transcripts for ncRNA is possible to do something similar to enrichment analysis? How can analyse the functions of the ncRNA that we get from the study?

Thank you very much in advance

Hi @Mee it might be that a RNAseq expert understands your concerns precisely, however could you provide a bit more background for everyone else please? For example, what package are you performing the analysis with, any description of the samples?


text:
How do I create a VCF file of all known pathogenic mutations in a gene of interest?
I would like to create a list (probably .vcf format would be good) of all known pathogenic missense mutations in a human gene of interest and then add other variants that could lead to the same pathogenic amino acid substitutions. I intend to search targeted sequencing data of my gene of interest for variants that would cause these particular substitions.
I'm very new to bioinformatics so I'm not so familiar with the appropriate databases and tools to use, but this is a rough idea of how it could be done:

extract all mutations listed in Ensembl within my gene of interest that are labelled with HGMD_MUTATION  (is there an option to download a vcf file with all HGMD mutations within my gene? I've only managed to do so for individual mutations.)
annotate list with reference and mutant amino acids for every mutation in the list
add additional mutations (which may not be present in HGMD) which could lead to the same substitution. E.g. I'm interested in the substitution Proline (reference codon: CCG) to Leucine (TTR or CTN) at AA position 102, so I would like to list all hypothetical mutations that would create a leucine codon in this position, not only the one given in Ensembl (which is CCG -> CTG).

Could you point me towards the best tools to use for these purposes? Unfortunately I haven't got a HGMD subscription. Many thanks!

Why HGMD specifically (not clear to me)? You can simply use the ClinVar database to find all pathogenic mutation within regions of interest/gene (you can check the genomic coordinates of any gene using either UCSC or ensemble).

I'm not sure how to generate the additional mutations, but I would say that HGMD is not the way to find all the pathogenic variants. I would probably filter this table by either Clinical significance or Evidence->Phenotype association.


text:
Scatter tasks around predefined batches inside a WDL workflow
Suppose I have two arrays Array[File] X and Array[String] Y both of length N, with Y being a sort of R-like factor that defines batch membership for X, for a total of M batches with M < N. A valid example would be:
Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

In which case N=5 and M=2. I want to scatter M tasks so that each task takes only the files corresponding to the same batch. In the previous example, one task would take as input:
Array[File] X =["file1", "file2", "file4"]
String Y = "A"

And the other task would take as input:
Array[File] X =["file3", "file5"]
String Y = "B"

Is there an easy way to do this within a WDL workflow in the general case? I understand that there are less elegant ways to do this as its own task but, as this seems like it could be a general purpose scatter/gather problem, an easy solution might be of general interest to the community.

Here's an alternative to using a custom task, albeit still rather unappetizing:
workflow groupby {
  input {
    File? nil
  }

  Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
  Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

  scatter (y in ["A", "B"]) {
    scatter (p in zip(X,Y)) {
      File? x_or_nil = if (p.right == y) then p.left else nil
    }
    Pair[String,Array[File]] group = (y, select_all(x_or_nil))
  }

  output {
    Array[Pair[String,Array[File]]] ans = group
  }
}

(It'll be slightly cleaner with the upcoming None literal instead of dummy nil)
FWIW, here's a lengthy discussion about the idea of allowing developers to define custom library functions which just run a small task that manipulates data structures (hinting that it needn't necessarily offload to a worker node, with the huge latency that typically entails).

When you use y in ["A", "B"] you are doing something that assumes you already know the unique contents of Y. This could be achieved using the keys() function but that is also only present in the development version and if the development version is to be used the solution with collect_by_key() is more minimalist.

It turns out the development version has a function collect_by_key that is particularly helpful for scattering around predefined batches:
Array[String] X  =  ["A",     "A",     "B",     "A",     "B"]
Array[File] Y = ["file1", "file2", "file3", "file4", "file5"]
scatter (pair in as_pairs(collect_by_key(zip(X,Y)))) {
  call test {input: X = pair.left, Y = pair.right}
}

This will call task t first with input:
String X = "A"
Array[File] Y =["file1", "file2", "file4"]

And then with input:
String X = "B"
Array[File] Y =["file3", "file5"]

This is a readable solution to this problem allowed by collect_by_key. This function is not very well advertised (I had completely missed it) but kudos to who introduced it in the language. The only drawback is that collect_by_key is not available in version 1.0.


text:
Sophisticated algorithms in biological systems
I am looking for examples of (known) algorithms/data structures in biological systems. My current understanding is that these don't exist - or haven't been found - beyond a relatively low level of complexity, which would seem like strong evidence to me that evolution isn't able to come up with something as complicated (and multi-step) as, e.g., the Schönhage-Strassen algorithm.
I heard of Bloom filters for recognizing smells in fruit flies (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6304992/).

It would be helpful to know more about what you mean by "algorithms" or "data structures". Organismal morphogenesis could be looked at as an incredibly complex set of linked algorithms- people call it a "developmental program".

There are many examples of biological systems doing computation, according to my understanding:

Regulatory networks encode decision trees of significant complexity: wikipedia.
Alan Turing did significant work on organismal morphogenesis showing how simple chemical patterns of inducers can lead to complex outcomes, "Turing patterns".
All of organismal development, cell division, etc. are tightly regulated procedures of gene regulation and morphogenesis.
The inspiration for neural networks is brains.


text:
Comparing multiple treatments to multiple other treatments in edgeR for simple effects in a complex experimental design
I am working with a RNA-seq data set in maize that has a relatively complex design. There are two levels of treatment A (nitrogen fertilizer level in the field, high or low), two levels of treatment B (nitrogen nutrients in in vitro cultures, high and low) and two levels of treatment C (two time points of sampling), all with 3 reps.
> library(edgeR)
> load("KC_Raw.RData")
> y <- DGEList(counts = KCraw.data[,2:25])
> keep <- rowSums(cpm(y) > 10) >= 3
> targets <- data.frame(rownames=colnames(KCraw.data)[2:25] ,
+                       Time=rep(c(rep("2DIC",12),rep("5DIC",12))) ,
+                       FieldN=rep(c(rep("FH",6), rep("FL",6)),2) ,
+                       CultureN=rep(c(rep("CL",3),rep("CH",3)),4))
> Group <- factor(paste(targets$FieldN,targets$Time,targets$CultureN,sep="."))
> targets <- cbind(targets,Group=Group)
> targets
   rownames Time FieldN CultureN      Group
1   KC1_H2L 2DIC     FH       CL FH.2DIC.CL
2   KC2_H2L 2DIC     FH       CL FH.2DIC.CL
3   KC3_H2L 2DIC     FH       CL FH.2DIC.CL
4   KC4_H2H 2DIC     FH       CH FH.2DIC.CH
5   KC5_H2H 2DIC     FH       CH FH.2DIC.CH
6   KC6_H2H 2DIC     FH       CH FH.2DIC.CH
7   KC7_L2L 2DIC     FL       CL FL.2DIC.CL
8   KC8_L2L 2DIC     FL       CL FL.2DIC.CL
9   KC9_L2L 2DIC     FL       CL FL.2DIC.CL
10 KC10_L2H 2DIC     FL       CH FL.2DIC.CH
11 KC11_L2H 2DIC     FL       CH FL.2DIC.CH
12 KC12_L2H 2DIC     FL       CH FL.2DIC.CH
13 KC13_H5L 5DIC     FH       CL FH.5DIC.CL
14 KC14_H5L 5DIC     FH       CL FH.5DIC.CL
15 KC15_H5L 5DIC     FH       CL FH.5DIC.CL
16 KC16_H5H 5DIC     FH       CH FH.5DIC.CH
17 KC17_H5H 5DIC     FH       CH FH.5DIC.CH
18 KC18_H5H 5DIC     FH       CH FH.5DIC.CH
19 KC19_L5L 5DIC     FL       CL FL.5DIC.CL
20 KC20_L5L 5DIC     FL       CL FL.5DIC.CL
21 KC21_L5L 5DIC     FL       CL FL.5DIC.CL
22 KC22_L5H 5DIC     FL       CH FL.5DIC.CH
23 KC23_L5H 5DIC     FL       CH FL.5DIC.CH
24 KC24_L5H 5DIC     FL       CH FL.5DIC.CH

I have used edgeR in R to calculate differential expression for contrasts involving 3 reps at one treatment combination to 3 reps at another treatment combination, for example
> y <- DGEList(counts = KCraw.data[keep,2:25], group = Group)
> y <- calcNormFactors(y)
>
> TMM <- KCraw.data[keep,2:25]
> for (i in 1:24) {
+   TMM[,i] <- TMM[,i] / (y$samples$lib.size[i] * y$samples$norm.factors[i]) * 1e6
+ }
>
> y <- DGEList(counts = TMM,group = Group)
>
> design <- model.matrix(~0+Group)
> colnames(design) <- levels(Group)
> y <- calcNormFactors(y,method = "TMM")
> y <- estimateDisp(y,design)
> fitQL <- glmQLFit(y,design)
> fit <- glmFit(y,design)
> myKC.contrasts <- makeContrasts(
+   H2H.H2L = FH.2DIC.CH - FH.2DIC.CL,
+   L2H.L2L = FL.2DIC.CH - FL.2DIC.CL,
+   H2H.L2H = FH.2DIC.CH - FL.2DIC.CH,
+   H2L.L2L = FH.2DIC.CL - FL.2DIC.CL,
+   H5H.H5L = FH.5DIC.CH - FH.5DIC.CL,
+   L5H.L5L = FL.5DIC.CH - FL.5DIC.CL,
+   H5H.L5H = FH.5DIC.CH - FL.5DIC.CH,
+   H5L.L5L = FH.5DIC.CL - FL.5DIC.CL,
+   H2H.L2L = FH.2DIC.CH - FL.2DIC.CL,
+   H5H.L5L = FH.5DIC.CH - FL.5DIC.CL,
+   H5L.H2L = FH.5DIC.CL - FH.2DIC.CL,
+   H5H.H2H = FH.5DIC.CH - FH.2DIC.CH,
+   L5L.L2L = FL.5DIC.CL - FL.2DIC.CL,
+   L5H.L2H = FL.5DIC.CH - FL.2DIC.CH,
+   levels=design)
> design
   FH.2DIC.CH FH.2DIC.CL FH.5DIC.CH FH.5DIC.CL FL.2DIC.CH FL.2DIC.CL FL.5DIC.CH FL.5DIC.CL
1           0          1          0          0          0          0          0          0
2           0          1          0          0          0          0          0          0
3           0          1          0          0          0          0          0          0
4           1          0          0          0          0          0          0          0
5           1          0          0          0          0          0          0          0
6           1          0          0          0          0          0          0          0
7           0          0          0          0          0          1          0          0
8           0          0          0          0          0          1          0          0
9           0          0          0          0          0          1          0          0
10          0          0          0          0          1          0          0          0
11          0          0          0          0          1          0          0          0
12          0          0          0          0          1          0          0          0
13          0          0          0          1          0          0          0          0
14          0          0          0          1          0          0          0          0
15          0          0          0          1          0          0          0          0
16          0          0          1          0          0          0          0          0
17          0          0          1          0          0          0          0          0
18          0          0          1          0          0          0          0          0
19          0          0          0          0          0          0          0          1
20          0          0          0          0          0          0          0          1
21          0          0          0          0          0          0          0          1
22          0          0          0          0          0          0          1          0
23          0          0          0          0          0          0          1          0
24          0          0          0          0          0          0          1          0
attr(,"assign")
[1] 1 1 1 1 1 1 1 1
attr(,"contrasts")
attr(,"contrasts")$Group
[1] "contr.treatment"

> myKC.contrasts
            Contrasts
Levels       H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
  FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
  FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
  FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
  FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
  FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
  FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
  FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
  FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
            Contrasts
Levels       L5H.L2H
  FH.2DIC.CH       0
  FH.2DIC.CL       0
  FH.5DIC.CH       0
  FH.5DIC.CL       0
  FL.2DIC.CH      -1
  FL.2DIC.CL       0
  FL.5DIC.CH       1
  FL.5DIC.CL       0

After analyzing these contrasts, I wanted to estimate some sort of simple effect, such as the culture media nitrogen level. To do this, I ran the following code.
> myKC.contrasts <- cbind(myKC.contrasts,
+                         Development = c(1,1,-1,-1,1,1,-1,-1),
+                         FieldN = c(1,1,1,1,-1,-1,-1,-1),
+                         CultureN = c(1,-1,1,-1,1,-1,1,-1)
+ )
> myKC.contrasts
           H2H.H2L L2H.L2L H2H.L2H H2L.L2L H5H.H5L L5H.L5L H5H.L5H H5L.L5L H2H.L2L H5H.L5L H5L.H2L H5H.H2H L5L.L2L
FH.2DIC.CH       1       0       1       0       0       0       0       0       1       0       0      -1       0
FH.2DIC.CL      -1       0       0       1       0       0       0       0       0       0      -1       0       0
FH.5DIC.CH       0       0       0       0       1       0       1       0       0       1       0       1       0
FH.5DIC.CL       0       0       0       0      -1       0       0       1       0       0       1       0       0
FL.2DIC.CH       0       1      -1       0       0       0       0       0       0       0       0       0       0
FL.2DIC.CL       0      -1       0      -1       0       0       0       0      -1       0       0       0      -1
FL.5DIC.CH       0       0       0       0       0       1      -1       0       0       0       0       0       0
FL.5DIC.CL       0       0       0       0       0      -1       0      -1       0      -1       0       0       1
           L5H.L2H Development FieldN CultureN
FH.2DIC.CH       0           1      1        1
FH.2DIC.CL       0           1      1       -1
FH.5DIC.CH       0          -1      1        1
FH.5DIC.CL       0          -1      1       -1
FL.2DIC.CH      -1           1     -1        1
FL.2DIC.CL       0           1     -1       -1
FL.5DIC.CH       1          -1     -1        1
FL.5DIC.CL       0          -1     -1       -1

Once I rerun the analysis for the CultureN contrast and look at the result for a particular gene, I see that it's estimated log2FC is equal to the sum of every simple contrast.
> lrt <- glmQLFTest(fitQL, contrast=myKC.contrasts[,"CultureN"])
> topTags(lrt,n=nrow(y$counts))["GRMZM2G445575",]
Coefficient:  1*FH.2DIC.CH -1*FH.2DIC.CL 1*FH.5DIC.CH -1*FH.5DIC.CL 1*FL.2DIC.CH -1*FL.2DIC.CL 1*FL.5DIC.CH -1*FL.5DIC.CL
                 logFC   logCPM        F       PValue          FDR
GRMZM2G445575 -6.63617 5.417106 151.5261 3.691525e-11 2.825777e-08
# FC is a data frame of the logFC of each constrast in columns for each gene in rows
> sum(FC["GRMZM2G445575",c("H2H.H2L","L2H.L2L","H5H.H5L","L5H.L5L")])
[1] -6.636197

My first question is if this analysis is a valid way of summarizing the simple effects of each treatment.
I would like to be able to also include the effects of the H2H.L2L and H5H.L5L contrast in the FieldN and CultureN comparison, but I am not sure how to do this, or if this would be valid because each of these contrasts includes treatments that have different levels of two treatment factors.

I don't quite get why you need to manually change the counts with a for loop. It is used as a size factor in the linear model. Now, for the pairwise comparison, what you have in the first part is ok

For the general culture N etc, run a linear model and run the estimateGLM... etc stuff, it should give you the effects you want. You cannot get everything from one model

I think the problem is in the design. There is no room for error or variation from a common ground.
I'm not sure of the solution, but I think you'll need to delete a column to give the design a degree of freedom, and probably add the intercept. You'll need to adjust the contrasts accordingly.
However, I would recommend to ask this at support.bioconductor.org. There are more experts on linear modelling and contrasts there than me. (If you ask there make it easier to copy and paste your code)

Thanks for the replies. I'm not really sure how to do the modifications you recommended to my design or contrasts. I am writing a question on bioconductor, but how can I make it easier to copy and paste my code? This was my first post on a bioinformatics support forum.

Ideally you could use a package as reprex. But if you use code blocks and remove the > at the beginning of the code, comment the output (so that if I paste it on the console it doesn't error) you'll do good. Another thing to consider is that I cannot load("KC_Raw.RData") because I don't have it, so try to make the example self contained. But the code provided is enough to understand your problem so don't worry if you use the same code used here there.


text:
.SDF file to .PDB file
I am currently working on identifying effective antiviral compounds against the corona virus. I have downloaded the Covid-19 AntiViral Candidate Compounds Dataset from the www.cas.org webpage, the link to which I have provided below. Me and a few other doctors are using Docking software to predict the docking capacity of these anti-viral compounds. Docking software requires the anti-viral compounds to be in the .pdb file format only then is it recognized by Automated docking software. Is there any way I can convert the .sdf file containing the anti-viral compounds to .pdb file or is there a way of obtaining the .pdb files for list of anti-viral compounds.
Thank you for your time and patience,
Anas Hussain.
Anti-Viral Candidate Compounds Dataset
https://www.cas.org/covid-19-antiviral-compounds-dataset

Use OpenBabel.
Something like babel -isdf file_in.sdf -opdb file_out.pdb.

Bizarrely, Obabel seems to have written a PDB that was incorrectly spaced for the OP as seen in Error pasing the following line in pdb

If you can use python, rdkit is a great package for cheminformatics —best installed via conda though.
It allows you do all sorts of things, such as match substructures, embed conformers, edit compounds, compute advanced parameters and many forms of shape & colour finger-/foot-printings.
Given a multi-compound sdf you can read each individual molecule and write them out as PDBs.
from rdkit import Chem
with i, mol in enumerate(Chem.SDMolSupplier('in.sdf')):
    Chem.MolToPDBFile(f'out_{i}.pdb')

We can be fancier with the output name
from rdkit import Chem
import re
with mol in Chem.SDMolSupplier('in.sdf'):
    # The headerline of a mol file is retrieved oddly.
    name = mol.GetProp('_Name')
    # make filename safe
    slugify = lambda name: re.sub(r'[\W_.-]+', '-', name)
    Chem.MolToPDBFile(f'{slugify(name)}.pdb')

To specify what name (e.g. CA) each atom gets and what is the residue name, chain and residue position —information absent in the sdf, it gets a bit odd.
  names = [' CA ', ' CB ', ' CG ', ' CD '] # say
  for i, atom in enumerate(mol.GetAtoms()):
        info = Chem.AtomPDBResidueInfo(names[i], # or whatever.
                                       residueName='LIG')

For more see this blog post.

for the rdkit sdf to pdb conversion, should i enter the multi entry sdf filename in the line below;
"with mol in Chem.SDMolSupplier('in.sdf'):"

Correct in.sdf is an example filename —the quotes must be there or it will think it a variable. If you are unfamiliar with python this route is probably not the best way.

I am unfamiliar with python however I am currently trying to learn it from scratch.
Thank you