text:
Is the quality score of fastq used somewhere besides trimming/fastqc?
In the fastq format every 4k+4-th line contains the positionwise qualityscore (ascii encoded):
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
...

which program uses this information besides trimming tools (cutadapt) or quality control (fastqc)?
I heard that segemehl (mapping) abandons this quality string in favor of running time (maybe other mapping tools too?)
But what about abyss, trinity, ...?
UPDATE:
So i atrifically generated a fastq file containing 2 identical reads that perfectly match a sequence in some genome (C.fna).
The two reads have the most extreme phred scores (only '!!'=highest sanger phread and only 'II'=lowest sanger phred)
Now using segemehl the mapping quality is identical (60)
@HD VN:1.0
@SQ SN:C_fna    LN:15440
@RG ID:A1   SM:sample1  LB:library1 PU:unit1    PL:illumina
@PG ID:segemehl VN:0.3.4    CL:segemehl -i C2.fna.idx -d C2.fna -q test.fq
A   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
B   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
[SEGEMEHL] Thu Feb  6 08:35:15 2020: matching w/ suffixarray has taken 0.000000 seconds.

So segemehl does NOT use the read quality phred (basically only using the read sequence for mapping)

Assemblers will generally ignore base quality scores since it's assumed that the input has been heavily QCed already.

I think almost all alignment programs would use base quality as a factor to calculate mapping quality (See here). These would also be used in Variant Calling to pick reads where high quality bases confidently align as either the reference base or one of the alternate bases.

see update above: segemehl does not use phred

@Paul Some aligners use phred (e.g., bowtie2), others largely ignore it. Phred scores have changed a bit over time and these days there are only 5 possible values, so their utility has decreased as sequencers have gotten more accurate.

But it is kinda weird that in this benchmark (https://www.ecseq.com/support/benchmark.html) segemehl is more accurate than bowtie2,bwa,... but slower and more memory hungry. If bowtie2 uses phred, why does it not gain a significant advantage over segemehl?

LoFreq is a variant calling tool that specifically leverages per-base quality scores.
For long read technologies (PacBio SMRT, ONT), I think quality score usage is mostly limited to qc steps as well, such as filtering out poor reads either by mean quality score of the whole read (as in ONT's basecalling softwares and Filtlong) or by the lowest quality score within a window in a read (Filtlong).  I think most long read mappers ignore these values (for example minimap2 and ngmlr do and can take fasta for read inputs).


text:
Differential gene expression bias due to effect of an individual sample
I am analysing a human single cell RNA seq experiment, where we have 4 groups, four samples each. Data has been analysed using Seurat, with the canonical workflow. I have tried DE using various methods (Wilcox, MAST, LR), and for each of them the top genes (biggest fold change, lowest adjusted p) are genes over expressed by specific samples, e.g.:

Looking at the average expression for these genes per each subject, and plotting on a Boxplot, it does not seem to be a genuine difference (e.g. for the gene above):

Besides the subject, there is no other parameter (e.g. age, sex) that is biased across groups and that shows up as differentially represented in the PCA plot.
Is there a DE method that takes into account this behaviour and corrects for it?
How can these results be presented?
Besides increasing number of subjects (not really feasible), is there a solution to analyse these data?
Thanks a lot!

How did you scale the data (ScaleData or SCTransform)? You can usually regress out any source of variation in your data using vars.to.regress = "Group". In your case, "Group" would be the column of object@metadata containing the Group information.

Thanks for the reply, I did not do SCTransform, and with ScaleData I did not regress any variable. I do not think that scaling by "Group" will be useful: 1) I do not want to mask genuine changes between groups, 2) the DE methods use normalised data, not scaled data.

Regressing out group would eliminate most if not all differential expression between groups. Regressing out sample would be a slightly more reasonable approach but would likely end up being pretty similar

This paper might be interesting to you https://www.biorxiv.org/content/10.1101/2020.01.15.906248v1

Thanks for the paper @alanocallaghan!

Leave-one-out cross validation is an iterative method that leaves out one sample until each sample has been left out once. You can then identify the intersection of the differentially expressed genes from the sample space of repeated analyses to produce a set of genes that are consistently differentially expressed independent of individual sample biases.
more info here:
http://efavdb.com/leave-one-out-cross-validation/

Thanks a lot Drew J., this seems a valid approach! However, I do wonder why there is not a DE method addressing this issue directly. I am not a statistician, so I may have missed some theory behind this, though.

The statistical models generally used for DE involve population level assumptions that are applied to incredibly small sample sizes. The main problem is that sequencing is still prohibitively expensive for small labs to produce data from experiments that include sufficient numbers to address effect size and power.

@Charles if you found the answer useful, please up-vote it. I am trying to build up my reputation points on this forum specifically to enable voting. Thank you.

One potential problem here is that, as far as I can see, you're not taking the hierarchical structure of your scRNA-seq data into account. I.e. the fact that cells from the same patient / sample will be more similar than cells from different samples, even within the same group. This is also the problem addressed in the paper mentioned by @alan. Also, it usually makes more sense to do this kind of group vs. group comparison for each subpopulation (cluster/cell type) separately (assuming you have multiple subpopulations in your data and did some kind of clustering).
One way of doing this is by using a pseudobulk-level approach, where you basically aggregate (by either averaging or summing) your data over each sample, for each cluster separately. So you end up with one counts matrix per subpopulation with one row for each gene and one column for each sample. The sample columns are then the sums (or averages) of all cells belonging to that cluster-sample combination.
You can then use bulk-RNAseq methods such as edgeR or DESeq2 on the aggregated data to test differences between groups of samples for each cluster. The muscat package implements all these steps and has quite extensive documentation. So you could check that and test it out on your data.
The paper describing the method of muscat can be found here: https://www.biorxiv.org/content/10.1101/713412v2

source: Stack Exchange

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

sort: votes


text:
Is it possible for coronavirus or SARS to be synthetic?
I have heard several conspiracy theories regarding the origin of the new coronavirus, 2019-nCov. For example that the virus and/or SARS were produced in a laboratory or were some variant of Middle Eastern respiratory syndrome (MERS), shipped via laboratory workers.
I am well aware bioinformatics has debunked many conspiracy theories involving infectious diseases, an important one being polio vaccination programs in Africa were the origin of HIV, for example here. Likewise, numerous conspiracy theories on the man-made origins of HIV were similarly overturned using bioinformatic studies.
Do you have bioinformatics evidence that would debunk the current conspiracists about coronaviruses?

More relevant at: https://skeptics.stackexchange.com/

Thanks for editing the question @terdon. There was a minor publication raising the 'spooky' flag in this regard (which got them publicity but was irresponsible) and yes it can be directly analytically debunked if the hypothesis concerns sufficient amino acids (in this case it was). It is a significant amount of work because you need to generate a null distribution - i.e. everything  - and that's a lot of work. Secondly there is a concern that a minor publication just gets further publicity.

I read your answer back when you first posted it and it is a really great resource! probably more useful to those of us in the field than laypeople, but great neverheless!

The scenarios are impossible and would be laughable if they were not so serious. The evidence is in the phylogenetic trees. Its a bit like a crime scene when the forensics team investigate. We've done enough crime-scenes often going to the site, collecting the pathogen, sequencing and then analysis - (usually neglected diseases) without any associated conspiracy theories.
The key technical issue is coronaviruses are zoonoses , pathogens spread to humans from animal reservoirs and phylogenetic tree really helps understand the  how the virus is transmitted.
Trees

The key thing about all the trees are bats. Bats lineages are present at every single point of the betacoronavirus phylogeny (tree), both as paraphyletic and monophyletic lineages, one example is this tree of betacoronaviruses here. Meaning the nodes connecting the branches of the tree to the "master-branch", represent common ancestors and these were almost certainly bat-borne cornaviruses. This is especially true for SARS and - here bat-viruses are EVERYWHERE.
The tree here also shows that SARS arose on independently on two occassions, again surrounded by bat lineages and 2019-nCov has emerged separately at least once, again associated with bats.
Finally the tree below is a figure from BioRxiv Zhou et al (2020) "Discovery of a novel coronavirus associated with the recent pneumonia outbreak in humans and its potential bat origin" shows the 2019-nCov lineage is a direct descendent of a very closely virus isolated from a bat (RaTG13*). This is a really conclusive finding BTW.

Note, I don't normally present inline images, but it is such a nice finding (hint to reviewers) and BioRxiv is open access.
Conspiracy theory 1: laboratory made virus
Literally it would require someone passaging a new virus, with unknown human pathogenicity, and independently introducing all the earlier passages enmass across bat populations of China. They would then hope each lineage becomes an indepedent virus population before, then introducing the virus to humans. Thus when field teams of scientists go around using mist nets to trap the bats, buy them from markets, isolate the virus and sequence it they would find a beautiful, array of natural variation in the bat populations leading up to the human epidemics, that perfectly matches vast numbers of other viral zoonoses. Moreover, this would have to have happen substancially prior to SARS and 2019-nCov, because the bat betacoronaviruses have been known about prior both epidemics, viz. its simply not feasible.
Biological explanation
General Bats are a reservoir host to vast numbers of pathogens, particularly viruses, including many alphaviruses, flaviviruses, rabies virus and beleived to be important in ebolavirus (I don't know about this) and even important to several eukaryotic parasites. It makes sense, they are mammals, so evolutionary much closer to us than birds for example, with large dispersal potential and roost in 'overcrowded' areas enable rapid transmission between bats.
Technical
The trees show bats are the common ancestor of betacoronaviruses in particular for the lineage leading into the emergence of 2019-nCov and SARS, this is seen in this tree, this one and the tree above. The obvious explanation is the virus circulates endemically in bats and has jumped into humans. For SARS the intermediate host, or possible "vector" was civet cats.
The theory and the observations fit into a seamless biological answer.
Conspiracy theory 2: Middle Eastern connection
I heard a very weird conspiracy theory attempting to connect MERS with 2019-nCov. The theory was elaborate and I don't think it is productive to describe here.
Biological explanation
All the trees of betacoronaviruses show MERS was one of the earliest viruses to diverge and is very distant from 2019-nCov, to the extent the theory is completely implausible. The homology between these viruses is 50%, so its either MERS or 2019-nCov. Its more extreme than mixing up yellow fever virus (mortality 40-80%) with West Nile virus (mortality <<0.1%), the two viruses are completely different at every level.
What about errors? Phylogeneticists can spot it a mile off. There are tell-tale phylogenetic signatures we pick up, but also we do this to assess 'rare' genetic phenomina. There is nothing 'rare' about the coronaviruses. The only anomaly is variation in the poly-A tail and that is the natural variation from in vitro time-series experiments. Basically we've looked at enough virses/parasites through trees, that have no conspiracy theories at all (often neglected diseases), and understand how natural variation operates - so a phylogenecist can shift the wheat from the chaff without really thinking about it.
Opinion
The conspiracy theories are deeply misplaced, and the only connection I can imagine is its China. However, the Chinese have loads of viruses, influenza in particular which causes major pandemics, but that is a consequence of their natural ecology (small-holder farming) allowing the virus to move between reservoir hosts. I've not visited small-hold farms in China, but I have in other parts of the world and when you see them, you get it. The pigs, chickens (ducks .. China), dogs, horses and humans all living within 10 meters of each other.
Conclusion
Shipping large numbers of bats to market, bat soup, raw meat from arboreal (tree-living) mammals such as civets that are sympatric to bats. Then consider the classical epidemiology in light of the phylogenetic data, which is very consistent, a single picture emerges that coronavirus is one of many zoonoses which has managed to transmit between patients.
Summary The fundamental point is the bioinformatics fit into the classical epidemiology of a zoonose.
*, Note The bat coronavirus RaTG13 predates the 2019-nCov outbreak by 7 years. It is not even clear whether the virus has been isolated, i.e. could just be a RNA sequence.
"They have found some 500 novel coronaviruses, about 50 of which fall relatively close to the SARS virus on the family tree, including RaTG13—it was fished out of a bat fecal sample they collected in 2013 from a cave in Moglang in Yunnan province."
Cohen, "Mining coronavirus genomes for clues to the outbreak’s origins" Feb, 2020 Science Magazine,

What about the "amino acid residues" in "4 inserts" that "have identity or similarity to those in the HIV-1 gp120 or HIV-1 Gag", that are mentioned by this Indian team, and which are probably the same mutations refered to by this Taiwanese professor? Did they make this up? Are their conclusions wrong?

We examined a very similar question here, https://bioinformatics.stackexchange.com/questions/11283/a-new-paper-suggests-the-corona-virus-has-uncanny-similarity-of-unique-inserts/11289#11289

Interestingly, a more recent paper suggests that the overall similarity with RaTG13 may be a little misleading, pointing to GD Pangolin-CoV instead: "Intriguingly, these same six critical AAs are identical between GD Pangolin-CoV and SARS-CoV-2 [16]. In contrast, although the
genomes of SARS-CoV-2 and RaTG13 are more similar overall, only one out of the six
functional sites are identical between the two viruses (Fig. 1B). [...]" https://doi.org/10.1093/nsr/nwaa036

That in turn seems to have caused controversy http://virological.org/t/response-to-on-the-origin-and-continuing-evolution-of-sars-cov-2/418

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


text:
Ensembl protein identifiers from different assemblies
In order to obtain orthologs for human genes, I am working with InParanoid version 7.0. This version of the database uses Ensembl Protein Identifiers (ENSP) from ENSEMBL version v54, which is based on NCBI 36 assembly of the human genome.
For my project, I need to have the Ensembl Gene Identifiers (ENSG) of these proteins, but in the last version of ENSEMBL, which is based on GRCh38 assembly. I need to use this version of InParanoid for consistency with previous analyzes.
Thus, my question is, could I just use biomaRt to map ENSP (from NCBI36) to ENSG (actual version) or this approach is inherently wrong?
Is it possible to "convert" ENSP identifiers from NCBI36 to GRCh38 and then do the map to ENSG?. I am aware of liftOver but I am not sure it would work here.

Those IDs are elderly! Ensembl 54 was 2009!
I would recommend using BioMart combined with the ID history converter. The ID history converter will convert old IDs to new, and BioMart will convert ENSPs to ENSGs. You can either use the ID history converter with the ENSPs first, then the current BioMart to get the ENSGs. Or you can use Ensembl 54 BioMart to convert the old ENSPs to old ENSGs, then update the ENSGs with the ID history converter. I would probably do both and compare results.

Totally agree! They are extremely old but I have to deal with them unfortunately. Thanks a lot for pointing me to the ID history converter, I wasn't aware about this tool. However, I think I didn't ask the question properly. Simplifying, could it be possible that a protein is ENSP00000387261 in v54 and then is ENSP00000387262 in v99? (Note that the example is hypothetical). I am not referring to the ENSPXXX.Y version of a protein but about if the identifier before the point could change between versions.

Yes, they could have new identifier numbers, which is why you need to combine the ID converter tool with BioMart. Most will stay the same, some will be dead and some will be different.

Perfect. I am trying with the tool. However, I am having troubles with some proteins. For instance, ENSP00000259791. Apparently it is retired since the release 75, according to its ENSEMBL entry. However, if you introduce it in the ID history converter and look at the table it is still matched to the same ID in the 99 release. If you download the results file, this is not happening anymore. Why this contradiction?. Am I confusing something or is it a bug? Thanks a lot for your help @Emily_Ensembl

That's weird. The tool should also give you three other IDs for that one, so you can go ahead with them. I will look into it.

Great. I am trying to parse the table for my 20,000 IDs in R so please let me know if you find what is the issue so I could take it into account systematically. Thanks!


text:
Seurat clusters
I'm trying to show cell cycle regression and input sample metadata to my umap.
I've been following Satija Lab's tutorials and have generated my umap but now want to specify which sample corresponds to each cluster, if one sample is driving a certain cluster, etc. I also, would like yo show cell cycle regression in each cluster. I following the cell cycle regression tutorial but that only tells me the board overview of the cells in each cycle. I would want to highlight it in the umap.
How would I go about each of these tasks?

You can use the group.by argument of the DimPlot() function. group.by accepts column names present in the meta.data column of the Seurat objects.

How do I add the sample and cell cycle regression information to the metadata?

Seurat::CellCycleScoring() stores the cell cycle information in the meta.data slot. Sample information can be added with Seurat::AddMetaData().

The other useful line for DimPlot (or FeaturePlot) is "split.by", which will make separate plots for each sample or whatever all side by side.


text:
Why does this naive Python paired Fastq reader work?
I wrote some throwaway code for ingesting paired-end reads from a Fastq file into memory.
def ingest(stream1, stream2):
    while True:
        record1, record2 = list(), list()
        for _ in range(4):
            record1.append(next(stream1))
            record2.append(next(stream2))
        yield record1, record2

I had initially put the for loop in a try/Except block, anticipating I'd need to explicitly break out of the endless while loop. However, the code seems to work correctly even after I removed it. I even "corrupted" the final read in one of the paired files, and the ingest function performs reasonably (the last read pair is never yielded).
My question is...why does this code work, since I never explicitly break out of the loop?
Here's an example of how the code is invoked.
>>> with open('seqR1.fastq', 'r') as fh1, open('seqR2.fastq', 'r') as fh2:
...     counter = 0
...     for r1, r2 in interleave(fh1, fh2):
...         counter += 1
>>> counter
5000
>>> r1
['@read5000/1\n', 'AGGTCTGGGGACAGGTGTTTTACAGGCTCCTTTAGGGGTGGCAAGTCTCCCGGAAACATGATTAGTGAAGACACGTTAACAGCGAAACTGACAGAGAACCTCAACATTCTCTGGTAGGAAAAGCCAGGATCCTGGCAGAGGAAGCAACGATTCGAGCAGGGAACTGTGCAGGACAGAGAACAATGCGGCCCCAATTTTGATCCAGGGGTAGTTAGCTTAGCTAAGGAGTTCCCCAGTCCCTAAGTGGAGAGGGGATCCAGAACCCTACAATTTATTATTATATCTTAGCGGATATATCTAT\n', '+\n', 'CCCCCGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG@GGGGGGGFAGGGCGGGGGFGGGGGGGGGCGGGGGDGGGGFGGGGGGGGGCGGGGGGGGGGGGGGGCGGGGGGGGGGGGGGGGG@G<GGGGGGGGGGGGFGGFGFF:DGGFFCGGGGCGGFEGGGGGGEGGGFB:GGGDG=FGGGG,GAGGGDGGFGGGCD3GFG+G?G*FGGG<GCGGGGC@CD<FFGF*:+G=FE<;GG,8C7GF=;9D=G<C5C66G6C>7)FD*6G7D)9;*C:)4D5+9>@G*=0):)*;1*504\n']
>>> r2
['@read5000/2\n', 'TGTGGGGTTCTGGATCCCCTCTCCACTTAGGGACTGGGGAACTCCTTAGCTAAGCTAATTACCCCTGGATCAAAATTGGGGCCGCATTGTTCTCTGTCCTGCACAGTTCCCTGCTCGAATCGTTGCTTCCTCTGCCAGGATCCTGGCTTTTCCTACCAGAGAATGTTGAGGTTCTCTGTCATTTTCGCTGCTAATGTGTCTTCACTAATCGTGTTTCCGGGAGACTTGCCACCCCTAAAGGAGCTTGTTGAACACTTGTCACCAGACCTAAGTCCTTAGGGATGAACCTGTGTCTCGGCAC\n', '+\n', 'CCCCBGC9FGCD-GGGFGGGEFEFF@CCGFGGGGFF+DG6,F69+FGG+FGF<FCC;GGEBFGGEDECFGGG6GDG7BFG9GCG9G9GDGGGG<GDFFFGEGC,<GC=GFG:,:A:EE+GD<:E,FGF=FCAGEFF++G,G,<EGGFGFFA>*FG9EGFCFCGC,B>43GFF@E@2>9ED@,568D=A;C+=5@E/F53D>,6F*9@+,E,C=*0+*@56<*9@:/F,1,:*)*14*0CC=*1)*0*,*)0;+)*))/**)**;0;//9)1+)10.()5/)070*)7*C**3/)(2*****\n']

Your code only runs prior to Python 3.7, now it raises RuntimeError: generator raised StopIteration. This is essentially a bug fix described in PEP 479 -- Change StopIteration handling inside generators. You can enable the new behaviour from Python 3.5 via from __future__ import generator_stop.
Note, it's not related to file-handles per se, you can reproduce it with a smaller example with a simple iterator:
it = iter([1,2])

def f(it):
    while True:
        yield next(it)

print(list(f(it)))

You should not write code relying on the bug

EDIT: Chris_Rands's answer is more authoritative and clarifies that this behavior is essentially a bug.  (I had only tested my example on Python 3.6.8 and hadn't seen PEP 479.)
When next runs out of lines, it raises StopIteration, and apparently that stops the while loop similarly to calling break without you needing to catch the exception.  I couldn't see a definitive source on that but you can see it in action if you catch and re-raise the exception:
def ingest2(stream1, stream2):
    while True:
        try:
            record1, record2 = list(), list()
            for _ in range(4):
                record1.append(next(stream1))
                record2.append(next(stream2))
            yield record1, record2
        except StopIteration:
            print("StopIteration raised")
            raise

You can also try to catch it outside the loop, but it doesn't get re-raised at that point.  (If you try raising one yourself in a while loop, though, it doesn't catch it in the same way, so there must be some extra Python magic going on.  I can't see any difference between the two exception objects so I'd be curious what's going on there.)  Neat example.  You could also make a pair of FASTQ parsers using Biopython's SeqIO.parse for more complex/permanent situations if needed.


text:
Can open-source software be peer-reviewed and published?
My colleague and I have developed a software tool intend to release it open-source. This tool is specifically for tasks in bioinformatics but we think it would be helpful for the wider community. Our institution will permit us to release it provided we get appropriate credit.
Thus we wish to publish it in peer-review. Is peer-review publication of bioinformatics software available? If so what is required to publish it?

There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck publishing in the Nucleic Acids Research yearly special webservers issue, since we make Galaxy wrappers around our tools.
You can also publish on tools in venues like Nature Methods, but it's vastly more difficult to get in there if your paper is purely on a tool and it's not something game-changing like salmon or kallisto.
Regarding what's actually needed to publish a tool, you will typically need the following:

Decent documentation. I ask for this as a reviewer and am asked about it when I submit things.
Easy installation. No one wants to spend two hours trying to compile your tool, so make sure there's a conda package or at least a docker container for it.
Test data. Without this your reviewers and users will likely not be happy since they won't be able to try it out on something small.

For the actual submission you'll tend to need comparisons to other tools. I'm generally not a big fan of this, but you will find that some reviewers will demand it (assuming there are other tools that do something similarish to yours).

Yes, there is even a journal dedicated to this: The Journal of Open Source Software
Here the software itself is peer reviewed and a decent documentation and tests are required. Since the submission and review process is based on GitHub, publishing there is straightforward if your code is on GitHub already.
Note: I am not affiliated with this journal, I just published two software packages there a while ago.


text:
Writing a python module to run Infernal on a list of fasta files
I am very new to using python3 and am hoping to use it to reduce the time I need to spend running Infernal to look for a variety of RNAs in a large collection of bacterial genomes. I am hoping to write a short module in python that will let me sequentially run infernal using a single covariance model on a list of genome files.
At the moment I am trying to use subprocess.run() or subprocess.Popen() to do so.

import subprocess
import os
import glob

genomes = list(glob.glob(os.path.join("/path/to/genomes/", "*.fa")))

for x in files:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", "genomes"], shell=False)

my issue is I do not know a way to have "genomes" passed as the components of the list within the command.
Is this a way off approach to trying to automate this? If so is there a better way to do it or is there a way to simply pass each file in the list to the cmsearch command?

what delimiter does cmsearch expect?  is the list separated with spaces, commas etc? 1 2 3 or 1,2,3 or 1, 2, 3 or something else?

The output of genomes = list(glob.glob(os.path.join("/path/to/genomes/", "*.fa"))) is comma delimted eg. /path/genome1.fa, /path/genome2.fa, ... For cmsearch it normal just takes a single file so  you would use the command: cmsearch /path/to/covariancemodel.cm /path/to/genomefile.fa            What I am hoping to do is replace that single genomefile.fa with each of the genome files in the genomes list. I hope that is clear! Thanks for the help.

You might be interested in a tool like snakemake.

I think what you want to do is:
import subprocess
import os
import glob

genomes = glob.glob(os.path.join("/path/to/genomes/", "*.fa"))

for gfile in genomes:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", gfile], shell=False)

each pass through the for loop will take a one of the genome files and pass it into the cmsearch command.
Note:
1) glob.glob returns a list so no need to cast it again - https://docs.python.org/2/library/glob.html

Thank you very much for the response! I ended up solving it a slightly different way with for x in genomes:  subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", "%s" %x], shell=False) but I just gave your recommendation and it works great and cleans things up a bit!

ok you shouldn't need "%s" %x since the values within the genomes list are already string.  Also str(variable) it the preferred way to convert and variable to a string


text:
How to calculate lipophilicity of proteins?
The question says it all. I want to calculate lipophilicity of proteins to screen the compounds out for antigen-antibody interactions.

Membrane affinity of protein is dependent on their structure and for a membrane associated protein only one side will be hydrophobic with a tell-tell charged ring around it (e.g. DHR1 domain). I.e. The structure or that of a homologue needs to be known.

I am not sure what you mean by " screen the compounds out for". Compounds normally means small molecules, but I have the feeling you mean the protein themselves (macromolecules). If so, that sounds like you are interested in doing global docking of a panel of protein vs. a threaded model (or similar) of an antibody, which is an absolute no go zone due to the propagated imprecision.

Agreed the question requires further clarification

A clarification:

I am not going to dock anything. I would like to calculate the lipophilicity of peptides, so still the question  'How?' remains. Also, protein is a chemical compound.

If it's oligopeptides see https://biology.stackexchange.com/questions/89177/what-is-the-role-of-the-pyrrolidine-piperidine-moieties-in-fagopyrins/89233#89233 for code

Oligopeptide
The lipophilicity of an oligopeptide or a small molecule can be calculated with RDKit, the classic compChem python toolkit. The measure is logP or partition coefficient, namely the magnitude of the ratio of the concentration in octanol over than in water, so the higher it is, the more it partitions in the octanol phase. However, with a compound of several dalton it may be good to see what the contribution of each atom to this score are. This measure is the Crippen contribution to logP, which can be handily plotted on the molecule. Note that these are may be at first glance look similar to Gastaiger partial changes, but are not.
Okay, first the imports
from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.Chem.Draw import SimilarityMaps

Now let's make the molecules. Let's start with something we can simply use a SMILES string of the web: Bortezomib.
name = 'Bortezomib'
mol = Chem.MolFromSmiles('O=C(N[C@H](C(=O)N[C@H](B(O)O)CC(C)C)Cc1ccccc1)c2nccnc2')

Now we can do all the relevant calculations on the mol object:
print('logP', Descriptors.MolLogP(mol))
AllChem.Compute2DCoords(mol)
contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, [x for x,y in contribs], colorMap='BuPu', contourLines=10)
fig.savefig(name+'_crippen.png', bbox_inches='tight')

LogP=0.3

colorMap='BuPu' is the matplotlib colorscheme. I have no idea which work best, I have spent way too much time trying and I do not know. But with BuPu, white is hydrophilic, purple is hydrophobic.
Say the oligopeptide does not have a SMILES or we don't want to do it that way. Let's take β-Casomorphins 1–3:
name = 'β-Casomorphins 1–3'
mol = Chem.MolFromFASTA('YPF')
mol = AllChem.ReplaceSubstructs(mol,
                                Chem.MolFromSmiles('N'),
                                Chem.MolFromSmiles('[NH][OH]'),
                                replacementConnectionPoint=0)[0]

The result with the last block from before gives:
logP=1

Proteins
With protein, it is a lot harder. Membrane affinity of protein is dependent on their structure. In PyMOL you can use the Adaptive Poisson-Boltzmann Solver (APBS) to show the surface by charge. For transmembrane protein there is clear ring of hydrophobicity —CHARMM-GUI membrane builder and a variety of tools for example can guess where it is. For a membrane associated protein only one side will be hydrophobic with a tell-tell charged ring around it (e.g. DHR1 domain) and there are no hard and tested methods for that.


text:
How to specify input resfile for Rosetta?
I'm trying to run the following Rosetta command.
rosetta_scripts.linuxgccrelease -s data/binding/dataset/test/protein.pdb \
    -parser:protocol data/binding/test/scripts/pack.xml \
    -out:path:all data/binding/test/output/ \
    -resfile data/binding/test/output/resfile.txt \
    -out:prefix RosettaOutput_ \
    -nstruct 1 -ignore_zero_occupancy false -overwrite

However, I keep getting the error:
core.pack.task.operation.TaskOperations: {0} [ ERROR ] File:resfile.txt not found!

ERROR: Cannot open file resfile.txt
ERROR:: Exit from: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc line: 1645
Error: {0} [ ERROR ] ERROR: Exception caught by JobDistributor while trying to get pose from job 'RosettaOutput_2MWY_0001'
Error: {0} [ ERROR ] Treating failure as bad input; canceling similar jobs

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/core/pack/task/operation/TaskOperations.cc:1645
[ ERROR ] UtilityExitException
ERROR: Cannot open file resfile.txt

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

protocols.jd2.FileSystemJobDistributor: {0} job failed, reporting bad input; other jobs of same input will be canceled: RosettaOutput_2MWY_0001
protocols.jd2.JobDistributor: {0} no more batches to process...
protocols.jd2.JobDistributor: {0} 1 jobs considered, 1 jobs attempted in 13 seconds

[ ERROR ]: Caught exception:

File: /home/sean/src/rosetta_src_2019.35.60890_bundle/main/source/src/protocols/jd2/JobDistributor.cc:327
1 jobs failed; check output for error messages

AN INTERNAL ERROR HAS OCCURED. PLEASE SEE THE CONTENTS OF ROSETTA_CRASH.log FOR DETAILS.

The error resolves itself if I move resfile.txt to the current working directory. Why can't resfile.txt be found?

A Rosetta question here? Excitement! Wait it’s solved. Darn.  For future questions... Rosetta commons forum is the ideal place for questions about Rosetta. Some users appear in both here and there, but there are more experts...

The problem is my pack.xml specified a resfile.txt location which over-wrote my command-line option.
    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile" filename="resfile.txt"/>
    </TASKOPERATIONS>

If you remove the filename option in the XML file, you are then able to specify the path.
    <TASKOPERATIONS>
       <ReadResfile name="Design_Resfile"/>
    </TASKOPERATIONS>


text:
Get nucleic acids' chain names from PDB API
There is a place in my pipeline where I call PDB's API with a ribosomal id to get all the subchains of that molecule programmatically
I'd like to separate subchains that are proteins from nucleic acids in the response, or at least know how to get the chain-ids of the nucleic acids so that i can exclude them from all manually.

Wondering if there is some query parameter like is_nucleic_acid in the PDB's ADL, but that's just a hunch. I'd be happy to learn if there is a more idiomatic way to do it( from Uniprot mappings maybe? ).
My query as of now looks like when i want to get, say uniprot acccessions for subchains(in fl):
const PDBSearchSelect = "https://ww.ebi.ac.uk/pdbe/search/select?"";
const params = {
// this is of the form "3j9m", "5afi", etc...
q: `${input.toLowerCase()}`,
//wondering if there is an "fl" parameter to get all nucleic acids
fl: "molecule_name, uniprot_id, uniprot_accession, struct_asym_id,...other",
//
rows: "200"
};

axios.get(PDBSearchSelect, {params}).... => ..

I would actually be super grateful if somebody could elucidate what's the difference between

https://www.ebi.ac.uk/pdbe/api/... and

`https://www.ebi.ac.uk/pdbe/search/pdb/select?

I've been looking at sides of the PDB API for a few weeks now and it still seems like complete babylon in terms of naming. I realize that it grew organically over years.
Thanks a ton.

Not sure I follow the first question you ask in your image about 'protein vs. NA' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/3j9m gives "molecule_type":"polyribonucleotide" and synonym":"16S ribosomal RNA" for "in_chains":["A"].

Also, the European PDB and the San Diego PDB have their own APIs each with different parameters. For some things one is better, for other the other. (or SIFTS). So sometimes its helpful to check what the other does. The San Diego PDB API is not as good and is documented in XML and in this case it is less good, but sometimes it is.

Namely whereas https://www.rcsb.org/pdb/rest/describeMol?structureId=3j9m has the RNA chain elements in this XML clearly have the attribute type="rna", the json (changing /rest/ to /json/ gives JSON, undocumented AFAIK) does not: https://www.rcsb.org/pdb/json/describeMol?structureId=3j9m.

@MatteoFerla, this is it, I think. You are the MVP.
Remember you from my other question, actually, and grateful to learn. The opaque documentation on many of these is quite frustrating at the same time.

@MatteoFerla, are there reasons that you know of other than historical as to why it all seems fragmented? Is that even true or am i just complaining?
I frequently dream of committing all of this to a huge intelligible GraphQL schema instead of dumping hours on trying to figure the right endpoint and query string every time i need to employ a new facet of a molecule.

No idea sorry. But I totally agree with your frustration. I needed to find the offset of the PDB residue numbers (not pose id) relative to Uniprot canonical transcript. SIFTS partially did it. I ended doing it manually, which was way quicker anyway...


text:
Calculating the mortality rate of a pandemic, e.g. coronavirus 2019-nCov?
When calculating how lethal a virus outbreak is, I've noticed that most sources use
CONFIRMED_DEATHS/CONFIRMED_INFECTIONS

However, given that there are only 2 possible outcomes; death or recovery, wouldn't it be more logical to use:
CONFIRMED_DEATHS/(CONFIRMED_DEATHS + CONFIRMED_RECOVERIES)

Which formula is applicable to which situations?

do you know that CONFIRMED_DEATHS CONFIRMED_INFECTIONS CONFIRMED_RECOVERIES are time function, and they will change for the time t

This may not strike most as a bioinformatics, but getting the key clinical outcome is essential in understanding the molecular basis of pathogenicity.
I think the mortality rate is over-reported. This is not to say the situation of 2019-nCov is not serious - it is very serious.
The two essential factors missing in your equation for 2019-nCov are:

Age. Mortality is likely skewed towards the old and very young
Seroprevalence. To identify the number of asymptomatic patients.

Old people have a weakened immune system (e.g. reduced T-cell response). The very young have an underdeveloped immune system. The problem is hospital cases are skewed towards those with severe infection, which is likely the old and very young. In 2019-nCov the first infecteds were middle aged men, but likely to be restricted to those exposed to the "Wuhan meat market".
Problem
Asymptomatic "spreaders" are known in influenza and are suspected in 2019-nCov. Asymptomatics will inflate the mortality rate (because we're not detecting them).
Solution(s)

One solution is to conduct seroprevalence surveys. Serological surveys screening for antibodies (IgG) against 2019-nCov will identify asymptomatic patients, using the first formula you stated. However, would you seriously go around the general population in a 2019-nCov epidemic taking blood samples? Errhhh no.
The bio-statistics approach is to use a catalytic model and critical to this is the age of the patient. Catalytic models are a priori models that incorporates the age-prevelance profile of infection. By introducing the incidence and age of hospitalisation of e.g. the old and the young, the rest of the distribution, i.e. asymptomatics can be calculated by regression. This estimates the total population of infecteds and again you use equation 1 to calculate the mortality rate recorded by the hospital.
Mathematical models based on transmission dynamics. The problem here is the rate of transmission is unknown.

I would assume the hospital is the place virus-induced mortality is reasonably accurate.
The second formula you described works very well for Ebolavirus, no-one has ever heard of asymptomatic Ebolavirus. It would also work well for another betacoronavirus Middle Eastern respiratory syndrome. The key problem with 2019-nCov is we don't know the number of asymptomatic patients or patients with such mild infections they never report illness.
Summary
In summary, I think the mortality is being over-reported if it is based on your equation 2. Equation 1, however, will also over-report the mortality rate.

i think that data is changing and no one can report the right mortality

To add to Michael G.'s answer, there is a further wrinkle in the case of 2019-nCov that is due to the exponential growth in confirmed infections daily.
That growth affects your second formula CONFIRMED_DEATHS/(CONFIRMED_DEATHS + CONFIRMED_RECOVERIES since deaths probably occur earlier on in the disease trajectory than recoveries do. Since the number of confirmed cases is growing exponentially on a daily basis, the number of deaths will presumably be inflated relative to the number of recoveries at this point.
This feature may have the opposite effect on your first formula CONFIRMED DEATHS/CONFIRMED INFECTIONS, depending on when the infections are confirmed. If the infections are generally confirmed several days before death occurs, then a large portion of confirmed infections are still too early to have ended in death and CONFIRMED DEATHS/CONFIRMED INFECTIONS actually underestimates the actual number of fatalities.

I was facing this problem in January when coming up with a methodology to forecast on the UrbanSurvival.com site.
My solution - using low data collection was to use a daily rolling rate of infections and divide that into the 5-day previous reported death counts.
This was based on the assumption that, for a given cohort, the people who get sick today will "get dead" in about 10-days time.
So half of this (5 days) ought to get me into the ballpark.
In Feb, when CDC was reporting 3.4% mort. rate, we'd been reporting that for a couple of weeks.
As a formula it's not unlike deriving an exponent.
That is:
Cases Today (T1) / Deaths 5-Days Ago (T2)
Unfortunately, this is presently running a bit higher mortality rate now (4.5%) but I would expect on the backside of the expansion of cases it will decline.
Here's hopin'...


text:
WGCNA: Problem with selecting soft threshold
I have 25 tumor samples with counts data. Initially, I filtered out low expressed genes and then converted counts to varianceStabilizingTransformation using deseq2 package. With this data I started using WGCNA for co-expression network analysis.
For selecting the soft threshold I see very strange plot. R2 cutoff is 0.8 and I see that none of the scale free topology model fit is above that.
Here is the code I used:
df is a dataframe with genes as rows and 25 samples as columns with counts data.
library("DESeq2")
filtered.counts <- df[rowSums(df==0)<5, ]

U3 <- as.matrix(filtered.counts)
vsd <- vst(U3, blind=FALSE)

oed <- vsd

gene.names=rownames(oed)
trans.oed=t(oed)
dim(trans.oed)

n=16462;
datExpr=trans.oed[,1:n]
dim(datExpr)

SubGeneNames=gene.names[1:n]

library(WGCNA)
options(stringsAsFactors = FALSE);
allowWGCNAThreads()

powers = c(c(1:10), seq(from = 12, to=20, by=2));
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
                      powerVector = powers,corFnc = cor,
                      corOptions = list(use = 'p'),networkType = "unsigned")

   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1 0.000273  0.041          0.786 4140.00  4000.000 6730.0
2      2 0.428000 -1.130          0.852 1570.00  1400.000 3770.0
3      3 0.673000 -1.540          0.882  729.00   579.000 2410.0
4      4 0.737000 -1.720          0.891  383.00   268.000 1670.0
5      5 0.745000 -1.830          0.886  220.00   134.000 1210.0
6      6 0.704000 -1.990          0.860  134.00    71.700  909.0
7      7 0.737000 -1.980          0.890   85.90    40.300  701.0
8      8 0.742000 -2.020          0.903   57.30    23.500  551.0
9      9 0.733000 -2.090          0.917   39.40    14.200  441.0
10    10 0.750000 -2.080          0.934   27.80     8.810  357.0
11    12 0.770000 -2.080          0.952   14.80     3.670  243.0
12    14 0.775000 -2.090          0.951    8.43     1.660  171.0
13    16 0.397000 -2.820          0.459    5.06     0.801  124.0
14    18 0.409000 -2.770          0.474    3.18     0.406   91.4
15    20 0.421000 -2.720          0.490    2.06     0.215   68.8

# Plot the results
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n", main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");

Is this the right way? Which soft threshold power should I select?
Any help is appreciated. thanq.

Don't worry about it too much. I would go with power 8 based on my general experience (also reflected in WGCNA FAQ) and on the mean connectivity around 50 and median around 20, which seems reasonable to me.

thanq very much Peter.


text:
Detection of CNV(InDel) of intermediate size
My impression is that small InDel (a couple of bp) is identified through cigar string in BAM and typical CNV (at least thousands of bp) is detected through read depth.
What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth.
Do we have consensus that this kind of CNV is hard to detect? Or is there a common strategy of detection for them? Do we generally think them having impact on protein function (especially in oncology background) ?

My impression is that small InDel (a couple of bp) is identified through cigar string in BAM...

Well, yes, by gaps in the alignment, easily detected by examining the CIGAR string.

...and typical CNV (at least thousands of bp) is detected through read depth.

Also by deviations from expected read pair distance or orientation.

What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth. Do we have consensus that this kind of CNV is hard to detect?

Yes, what you describe has certainly been recognized and acknowledged in the literature. It was one of the primary motivations for a recently published paper I contributed to on de novo (germline) variant discovery. Here's an excerpt from that paper:

In a reference-mapping context, calling indels with confidence requires accurate mapping of each read spanning the indel, with all gaps arranged consistently. This is possible only for short indels and tends to be prone to error and misalignment. Thus prediction of indels with length 10 bp has proved to be very challenging and accompanied by high false-positive and false-negative rates. Furthermore, the prediction of SVs via read mapping is only possible through indirect signatures such as alterations in read depth or read-pair signatures. These signatures can be quite noisy and result in high rate of false-negative and false-positive prediction.

As for your question:

is there a common strategy of detection for them?

Numerous mapping-free k-mer based methods for variant detection are coming up in the literature, and all seem to perform pretty well. The trick is how to define which k-mers are interesting—then you grab those k-mers (or the reads containing them) and assemble them into variant-spanning contigs. A few examples:

In the Kevlar paper, we defined  "interesting" k-mers as those with high abundance in an individual of interest but zero (or very low) abundance in the individual's parents.
The HAWK paper, which describes a k-mer based GWAS-like method, defines "significant" k-mers as those significantly over- or under-represented in one class of samples versus a second class.
The NovoBreak paper compares paired tumor and normal WGS data to find k-mers that span mutatino breakpoints in cancer genomes.

Do we generally think them having impact on protein function (especially in oncology background)?

There is definitely potential that indels can have a substantial impact on protein function. But since their detection has been so difficult for so long, a lot remains unknown about their actual collective impact on function.

There is a trade off between read length and the limit of the indel that can be accurately mapped. Intermediate-size indels specifically are better detected with longer read lengths.


text:
tp53 example in Momand's book
I am new to bioinformatics. While reading Momand's 'Concepts in Bioinformatics and Genomics' , I understood clearly that each codon translates to a specific amino acid.
But then, he talks about the TP53 gene, and then without explanation talks about P53... Is that a naming convention? It's not explicitly explained... Does that applies for any case? ie TP90 gene encodes protein P90?
I don't understand the naming: correct me if wrong, but they named the gene after the protein's weight, which is 53 kilodaltons?
(It confuses me with the location of the gene).
Lastly I do not understand: if there is a mutation in the DNA, the mutation shall appear in only one cell I guess, how then the others cells' DNA start to mutate? or is it that all DNA in all cells are mutated instantly?

This is possibly more of a molecular genetics question appropriate to biology.stackexchange.com, also note that it's usually better on stackexchange sites to restrict your post to a single question and do some preliminary research via google or your search engine of choice. But I don't see any harm in answering it here nonetheless.
The real answer is that naming conventions for genes/proteins are incredibly inconsistent and counterintuitive. You can't really draw any conclusions about a gene/protein based on its name alone. For example, some proteins such as p53 are named for biochemical features, others are named for phenotypes observed when their gene is mutated, such as for instance the gene decapentaplegic in fruit flies. Often they involve little jokes, such as the gene called mothers against decapentaplegic which represses the decapentaplegic mutant phenotype. So it's all a bit of a mess and people kind of assume that you'll follow along, which can be confusing.
It's a little rarer that gene names differ substantially from the names of the proteins they encode, but not unheard-of. This is indeed the case with your TP53/p53 example. See a short explainer here or a longer one here. What you'll notice is that there are in fact multiple different names for the same gene/protein. What you'll often see is that the shorter name (p53) is the one that people tend to use colloquially.
For your second question regarding how mutations are propagated, the short answer is that mutations are inherited from mother to daughter cells. Cells don't generally "infect" other cells with mutations (though in biology there are always exceptions). I think that what you are interested in is the difference between "somatic" and "germ-line" mutation.
Hope that helps.

My paperback edition is 'not for sale in the US' and lacks some clarifications on the margins which I have seen in an online version of the book. These clarifications are really good, not sure why they are missing!!

@Chicago1988 I think that these clarifications are missing from most textbooks in biology (!)


text:
How to translate amino acid sequences to Nucleotide sequences
I want to convert a list of fasta ( protein sequences) in a .text file into corresponding nucleotide sequences. A Google search gives me result of DNA to protein conversion but not vice versa. Also, I came across How do I find the nucleotide sequence of a protein using Biopython?, but this is what I am not looking for. Is there any possible way to do it using python.Moreover, I would like to solve it using python programming. I am sure there must be some way to do it rather than writing a code from scratch. Thanks!

Cross-posted: https://www.biostars.org/p/421074/

Although there is not a unique nucleotide sequence that translates to a given protein, one can list all the possible DNA sequences that do translate to that protein.
An online tool that does just that is Backtranambig, from EMBOSS. It produces a DNA sequence representing all the nucleotide sequences matching the input protein, using IUPAC ambiguity codes.

DNA Chisel (written in Python) can reverse translate a protein sequence:
import dnachisel
from dnachisel.biotools import reverse_translate

record = dnachisel.load_record("seq.fa")
reverse_translate(str(record.seq))
# GGTCATATTTTAAAAATGTTTCCT

You can't do this because there is redundancy in the genetic code and the same protein sequences can be encoded by different nucleotide sequences. There are 64 codons and ~20 amino acids, e.g. GCT, GCC, GCA and GCG all encode Alanine.

For the sake of the programming exercise, one could use IUPAC ambiguity codes to reflect the set of possible codons for a given amino-acid.

@bli Agreed, one could even build a model to predict the underlying codons since there is codon usage bias. For practical purposes I'd do a TBLASTN against NT or similar


text:
Error in createMAE function: non-unique values when setting 'row.names' in TCGA LIHC Data
I am working with the createMAE() function of the ELMER Bioconductor package. While executing the createMAE() function using TCGA LIHC DNA Methylation 450K and RNASeq HTSeq - Counts data, I am getting this error:
Error in.rowNamesDF <-(x, value = value) :

duplicate 'row.names' are not allowed

In addition: Warning message:

non-unique values when setting 'row.names': 'TCGA-DD-AACA-02A', 'TCGA-DD-AACA-02B'

Can anyone tell me how to get rid of such error?

The error message is pretty clear: R requires unique row names on data frames and some of the "values" that you are trying to set as row names are non-unique.
You will need to provide a vector i) of length the same as your row number, ii) that is composed of non-unique values.


text:
No output when trying to obtain protein sequence from PDB file
I am trying to obtain the sequence from PDB file, but no output is given when I run the following code :
from Bio import SeqIO

with open("1a9l.pdb", "r") as handle:
    for record in SeqIO.parse(handle, "pdb-atom"):
        print(record.seq)

If you only want protein sequences for an PDB code, NCBI has already a dataset (in genbank format) for that. Therein the accession IDs are the code underscore chain letter...

I was trying to do it through the script,as it is part of pipeline..

If you are using a conda version of python you can use Pymol module for this, which is actually more robust in some aspects than bio.PDB, although errors may cause segmentation faults. The following:
import pymol2
## paralellisable version.
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1a9l')
    print(pymol.cmd.get_fastastr('polymer.nucleic'))

Gives me:
>1a9l_A
GGGUGACUCCAGAGGUCGAGAGACCGGAGAUAUCACCC

Alternatively, you can use PDBe's API.
import requests

response = requests.get('https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/1a9l').json()['1a9l']

for entity in response:
    if entity['molecule_type'] == 'polyribonucleotide':
        print(entity['sequence'])

Do note that missing 'residues' (residues, small molecules and nucleotides) will cause issues. PDB files do not have missing 'residues', while mmCIF files have a dictionary that does. The API will give you the latter.

Hi, I tried conda install -c schrodinger pymol to install pymol module. On import, it shows me an error : ModuleNotFoundError: No module named 'pymol2'
Is this the way to download it for use with conda?

conda install -c schrodinger pymol will install both. Does import pymol work? That is, is the version of python you are running the conda one? If so, single threaded pymol in the example runs as import pymol without the context manager line. Doing conda search -c schrodinger pymol for me (Ubuntu) gives me choices including pymol 2.1.1 for  py27_2 and py36_2. Do you get the same versions?

yes,I get same.No import pymol is also not working,yes I am running on same version of python as conda.

Okay, just to clarify: (a) there is no sudo'ing/nohup/systemctl service going on, (b) your which python (or whichever you are calling) and your which conda point to the same folder and (c) when you type conda list pymol it gives it as installed, correct?

I should have phrased it not as which python, but within your run sys.executable gives that one.

yes,it shows:
Name                    Version                   Build  Channel
pymol                     2.3.4            py37h95b93ae_0    schrodinger

Installed with conda pymol is no different than numpy or pandas, so this is odd if the python is indeed correct one (cf. sys.executable ) and I have never had issues with it on dozens of systems (MacOS, Debian, RedHat) and architectures (x64, arm64) w/o issues. Uninstalling and reinstalling it?

Does pip freeze | grep 'pymol' - give the same result? Does print([k for k in sys.modules.keys() if 'pymol' in k.lower()]) gives you something? If so, pymol2 = importlib.import_module('pymol2') just might work, but is nearly semantically equivalent to import pymol2 so I highly doubt it.

I will try to troubleshoot this issue this week.If import issue still persists,will ask you.

That PDB entry is an RNA molecule, so there is no protein sequence to extract.
If it extracted the RNA sequence by default here you can see there could be an ambiguity as to whether it was an RNA sequence or a protein sequence, as the alphabets share characters.

Got it..So,I can't extract RNA sequence through this option.

Not using SeqIO from Biopython, no. The "pdb-seqres" option won't work here either. You could download the flat file of PDB sequences and read them into a Python dictionary: https://www.rcsb.org/pages/download_features#FASTA. Be aware though that sometimes atoms are missing in the ATOM records, which may have been why you wanted to read directly from the atoms in the first place.


text:
Why does the FASTA sequence for coronavirus look like DNA, not RNA?
I'm looking at a genome sequence for 2019-nCoV on NCBI. The FASTA sequence looks like this:
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
...
...
TTAATCAGTGTGTAACATTAGGGAGGACTTGAAAGAGCCACCACATTTTCACCGAGGCCACGCGGAGTAC
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

Coronavirus is an RNA virus, so I was expecting the sequence to consist of AUGC characters. But the letters here are ATGC, which looks like DNA!
I found a possible answer, that this is the sequence of a "complementary DNA". I read that

The term cDNA is also used, typically in a bioinformatics context, to refer to an mRNA transcript's sequence, expressed as DNA bases (GCAT) rather than RNA bases (GCAU).

However, I don't believe this theory that I'm looking at a cDNA. If this were true, the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail. But I believe the coronavirus has a poly(A) tail.
I also found that the start of all highlighted genes begin with the sequence ATG. This is the DNA equivalent of the RNA start codon AUG.
So, I believe what I'm looking at is the true mRNA, in 5'→3' direction, but with all U converted to T.
So, is this really what I'm looking at? Is this some formatting/representation issue? Or does 2019-nCoV really contain DNA, rather than RNA?

Excuse my ignorance, but I expected a lot more AAA (3' poly(A) tail) at the end of that sequence. Am I confusing multiple formats here or is your sequence missing a part? Is this the result of transcribing RNA as DNA?

should be 11 triplets, 33 times A

That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil).
The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse transcribed, either by targeted reverse transcription or random amplification and thus converted to cDNA. cDNA is stable and is essentially reverse transcribed RNA.
The cDNA is either sequenced directly or further amplified by PCR and then sequenced. Hence the sequence we observe is the cDNA rather than RNA, thus we observe thymine rather than uracil and that is how it is reported.

It is really hard to sequence directly from a patient It has been done (Direct RNA Sequencing of the Complete Influenza A Virus Genome) but of course that doesn't change the correct answer here - Virtually all sequencing of RNA viruses is done through cDNA.

Hi @iayork, thanks I accept its doable, but e.g. one risk is one sloppy wet-lab mistake and the vRNA is degraded, e.g. vRNA is stored at -70oC and needs RNase inhibitors, viral cDNA is stored at -20oC and no-one cares. Sometimes there are legal issues (these are pathogens).

Most sequencing experiments, be it Illumina-based next-generation-sequencing or Sanger sequencing uses DNA as template, not RNA. Even if this virus is RNA-based it would be reverse-transcribed prior to any sequencing experiment. Therefore the output is DNA and this is what NCBI provides here.

If this were [cDNA], the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail.

A cDNA sequence, maybe confusingly, refers to the coding strand of the cDNA (despite being called “complementary”). So while cDNA is the result of reverse transcribing RNA into DNA, by convention it has the same strandedness as the original RNA. That’s why what you’re seeing is read in 5′→3′ direction and contains a visible poly(A) tail. Having a single conventional reading direction for all archived sequences vastly simplifies data handling, and reduces errors.
In fact, since cDNA is double-stranded, there is no a priori reason why a computer-stored cDNA sequence should refer to the template strand (i.e. the opposite strand, which is synthesised from the RNA during reverse transcription).
The whole (simplified) synthesis process of cDNA is as follows:

A primer hybridises to the template RNA molecule.
The RNA template is reverse transcribed into DNA using reverse transcriptase.
The RNA template is removed.
A complementary strand is transcribed along the (currently) single-stranded cDNA, resulting in a double-stranded cDNA product.

It's not common to sequence directly from RNA because most sequencing platforms don't have that as an option. Nanopore sequencers do allow this, but I'm not aware yet of any 2019-nCov preprints involving nanopore RNA sequencing. I expect that will change in the next month or so.
Commercial kits exist; there are no insurmountable technical issues with it. Direct RNA sequencing can be done locally, on-site near the point of discovery without sample transfer or culture on a USB-powered device that fits in a pocket (RNA preparation takes about 2 hours). Flow cells that have potentially-infectious RNA inside them can be disposed as biohazard waste. However, the ease at which RNA can be quickly converted to more stable cDNA then amplified to create a much higher concentration DNA sample (which is quicker / more efficient to get results from) means that cDNA is generally preferred for sequencing unless the native RNA is needed (e.g. for looking at RNA base modifications that are destroyed when converting to cDNA).
There is a paper on coronavirus direct RNA sequencing with nanopore here; I would expect that 2019-nCoV would have a similar difficulty. The zika virus has an extremely low viral load in human blood, but has also been sequenced via direct RNA sequencing of [carefully] cultured cells (see here).
Regardless of whether or not RNA sequencing has actually been carried out, most genetic data analysis programs will only work with A/C/G/T sequences, so it's conventional to replace any U parts of an RNA sequence with T for data storage. There's no loss of information by doing this, as T replaces all Us in the RNA sequence.


text:
hisat2 --rna-strandness option and downstream htseq-count analysis
I've got some doubts on the hisat2 --rna-strandness option and its output for downstream analysis. Is it expected to see a difference in alignment and counting results given the default usage of --rna-strandness in hisat2 followed by htseq-count -s reverse on a strand-specific assay?
Please see below.
I understand that the --rna-strandness option produces an XS tag to indicate where a transcript is from (on the + or - strand) for downstream transcriptome assembly analysis. I have a paired-end stranded sequencing library that was aligned to the genome using hisat2 without specifying the --rna-strandness (in other words, the default unstranded was the usage). Following this, the reads were assigned to genes using htseq-count and this time -s reverse was specified given the strand-specific sequencing assay type.
Since --rna-strandness is for transcriptome assembly using the XS tags generated and htseq does not use XS tags for counting, I presume there should be no practical impact from the above. Could you also shed light on this? in case I may have been overlooking other facts of the usages of the tools.
To help verify the above, I re-aligned and counted the reads from 2 samples by switching on --rna-strandness RF in hisat2. I attach the alignment and count features info. below for assessment.
Overall alignment rate of Sample 1:
94.52% (--rna-strandness RF)
94.12% (--rna-strandness unstranded)

Overall alignment rate of Sample 2:
94.57% (--rna-strandness RF)
94.15% (--rna-strandness unstranded)

Feature counts of Sample 1 (following --rna-strandness RF + -s reverse):
__no_feature        6327294
__ambiguous     2954776
__too_low_aQual     3784481
__not_aligned       688856
__alignment_not_unique      4858182

Feature counts of Sample 1 (following --rna-strandness unstranded + -s reverse):
__no_feature        6291151
__ambiguous     2911298
__too_low_aQual     4075017
__not_aligned       754400
__alignment_not_unique      16136045

Feature counts of Sample 2 (following --rna-strandness RF + -s reverse):
__no_feature        5417882
__ambiguous     1708510
__too_low_aQual     3532352
__not_aligned       564596
__alignment_not_unique      2859501

Feature counts of Sample 2 (following --rna-strandness unstranded + -s reverse):
__no_feature        5359434
__ambiguous     1676091
__too_low_aQual     3813344
__not_aligned       623122
__alignment_not_unique      2891792

These results look comparable to me across pipelines.
Thanks
Guan

Maybe this post can help https://bioinformatics.stackexchange.com/questions/4074/hisat2-which-option-should-mention-for-strand-specific-library-read/4075#4075

If you reran the command with the correct settings, just leave it at that.  (It is not at all clear to me that strandedness rf is correct)
If you want people to tell you if you ran the commands right, you need to put down what commands you used.

Sample 2 looks comparable but Sample 1 shows a drastic reduction in unique alignments.
Your strandedness param is going to dictate how your alignments are done in the sense that discordant pairs will either be detected or not.
I would identify the strandedness of your library and use the appropriate respective param.


text:
Showing name of genes in boxplot
I have a long format data frame like this
> head(df)
     gene

1:  CCL28   2.030674       2.023333 1.921508    1.952007
2:  CXCL2   2.204847       2.275296 2.175790    2.241554
3:  CXCR4   2.121186       2.147181 1.950065    1.982639
4: IFNGR1   2.311390       2.324111 2.220727    2.255556

And by this code
ggplot(df, aes(variable,value)) +
  stat_boxplot(geom="errorbar", width=.5)+
  geom_boxplot(aes(fill=variable))+
  theme_bw()+
  theme(axis.title.x=element_blank(), axis.title.y=element_blank())+
  stat_summary(fun.y=median, colour="red", geom="line", aes(group = 1))+
  geom_jitter(position = position_jitter(0.2))

I have this box
I have added this line to my code where I got
geom_text(data = df_melt, aes(x = variable, y = value, label = gene), size = 5)

[![enter image description here][2]][2]
I want to show the name of genes in box plot but I don't how, cam you help?
Thanks

Surprisingly the first link on google for "add text to boxplot ggplot" gave the answer: https://stackoverflow.com/questions/23096108/ggplot2-add-text-on-top-of-boxplots

I think you need to edit your question because the colnames of head(df) and the df you use in ggplot2 are not the same (no mention of the column variable and value. So, I believed that you reshape your dataframe. Please indicate all step and provide a reproducible example

In addition, to have several errors in your code for plotting (for example, you are using df for plotting and then df_melt for geom_text but you didnot mention what is it). Based on your graph, it seems that you have issues with the melting part.
However, these are quite secondary issues that are probably be fixed with the code below. with your example, there is one trick because you are using jitter to plot your points (that is not addressed by the link that @ATpoint provided in comments). This will randomly position the point at the right y value but with different x position in order to not overlap points.
So, when you will add your text label, the labeling won't fit to the position of points. But you can make the jitter reproducible by adding a seed (see this github discussion: https://github.com/tidyverse/ggplot2/pull/1996). However, geom_jitter function does not have the argument seed, so you can use geom_point function and the argument position_jitter instead.
So, using the few genes that you provided as an example. I have to reshape it (a step you probably did but you did not provide in your question), here I'm using pivot_longer function from tidyr package:
library(dplyr)
library(tidyr)
DF <- df %>% pivot_longer(., -gene, names_to = "variable", values_to = "value")

# A tibble: 16 x 3
   gene   variable       value
   <chr>  <chr>          <dbl>
 1 CCL28  Res_immune      2.03
 2 CCL28  Non_res_immune  2.02
 3 CCL28  Res_OBP         1.92
 4 CCL28  Non_res_OBP     1.95
 5 CXCL2  Res_immune      2.20
 6 CXCL2  Non_res_immune  2.28
 7 CXCL2  Res_OBP         2.18
 8 CXCL2  Non_res_OBP     2.24
 9 CXCR4  Res_immune      2.12
10 CXCR4  Non_res_immune  2.15
11 CXCR4  Res_OBP         1.95
12 CXCR4  Non_res_OBP     1.98
13 IFNGR1 Res_immune      2.31
14 IFNGR1 Non_res_immune  2.32
15 IFNGR1 Res_OBP         2.22
16 IFNGR1 Non_res_OBP     2.26

And then, for the plotting part, you can use geom_text_repel function of the ggrepel package to add labeling without overlapping your points:
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

Warning
If you have a long list of genes, it will be pointless to show all of them on the boxplot. It will makes things really illisible. If you want to show some particular genes of interest, make a subset of the original dataframe. A lot of posts on SO and here have addressed this subsetting.
Reproducible data example
structure(list(gene = c("CCL28", "CXCL2", "CXCR4", "IFNGR1"),
    Res_immune = c(2.030674, 2.204847, 2.121186, 2.31139), Non_res_immune = c(2.023333,
    2.275296, 2.147181, 2.324111), Res_OBP = c(1.921508, 2.17579,
    1.950065, 2.220727), Non_res_OBP = c(1.952007, 2.241554,
    1.982639, 2.255556)), row.names = c(NA, -4L), class = c("data.table",
"data.frame"), .internal.selfref = <pointer: 0x563d91b65350>)

Did OP add comments and then delete them? This seems to be a one-way conversation and OP either deleted their comments or edited their posts, adding their comments in-post.

You're right. OP deleted his comments. It was mainly asking for explanation about error message that I presumed have been solved. I will remove my comments and they are not anymore relevant

That behavior is both entitled and ungrateful. Wow.


text:
Can I submit a R package to Bioconductor or CRAN if I have already published it a journal?
I have written a bioinformatics package in R that I want to publish in a bioinformatics Journal. Presently, I am maintaining a local repo of that package and I want to put in the Bioconductor repository (after publication or at least after submission to the journal), is it possible? or I must submit it to Bioconductor first?
Thanks

You can submit a package to a repo after the publication, but I would say that you really should do it before.
Bioconductor accept only packages that are not published on CRAN, however, academic publications are fine.
Conversely, Bioinformatics does not enforce any particular platform for sharing the code, but you have to make it available somehow (e.g. GitHub repo, university webpage, bioconductor). From the instructions for authors

If the manuscript describes new software tools or the implementation
  of novel algorithms the software must be freely available to
  non-commercial users at the time of submission, and appropriate test
  data should be made available.

However, you as a developer should make all the effort to make your software available and easily accessible. Also, sooner you make it available, sooner people start to cite you, which is a great advantage for you. Finally, you really need users to make sure that your package really works (i.e. debugging all the corner cases).
I could not emphasise more, the correct order is 1. make it available 2. write a publication

Thanks for your time. Actually, I wanted to whether I can still submit the package to Bioconductor even if its already published in a Bioinformatics journal.

Ah I see. I added also info from the Bioconductor side.


text:
how to get FlyBase ID conversion table
FlyBase uses systematic naming of genes, transcripts and proteins in a format FBtypeNUMBER. However, the numbers are not the same for genes and their products. For example gene FBgn0259101 transcribes to FBtr0299513 and translates to FBpp0288787.
On this webpage users can conver individual IDs, but what If I need to convert practically whole transcriptome? Where can I get the full translation table?

On the right hand side you can upload file with IDs. You can  upload a file with the whole transcriptome and convert to polypeptide, you can also upload from gene > transcript and then once morte to polypeptide and merge the 2 tables yourself later.

Indeed the web interface does not seem to link the full translation table, but it is available on their ftp server at releases -> release_ID -> precomputed_files -> genes -> fbgn_fbtr_fbpp_fb_2019_06.tsv.gz
For example, this is the latest (FB2019_06) translation table available.

took me like 20 minutes to find it, therefore I thought it might be worth sharing this with others.


text:
Griding ggplot object in r
I have this data frame

head(d)

# A tibble: 6 x 3
  gene  variable       value
1 CCT5  Res_immune     10.6
2 CCT5  Non_res_immune 10.5
3 CDK1  Res_immune      9.25
4 CXCL5 Non_res_immune  9.82
5 GBP1  Res_immune      9.04
6 GBP1  Non_res_immune  9.43
>

And this code
library(ggplot2)
library(ggrepel)
ggplot(d, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)

The genes in first column are involved in cytokine signaling pathway in two groups of patients
This is my plot

I have another data frame for genes in IL17 signaling in the sample patients like
head(d)
# A tibble: 6 x 3
  gene  variable       value
  <chr> <chr>          <dbl>
1 CDK6  Res_immune     17.6
2 IL6  Non_res_immune 14.5
3 CXCL8  Res_immune      1.25

How I can put plots from both data in a grade like below where the name of pathways shown instead of DEL, INS, SNP and total like below picture?

you should add a reproducible example of both dataframe (Cytokine and IL-17 signaling) by providing the output of dput(head(d)). It will make things easier for people trying to help you.

What you are looking for is facetting. Using this key word on any search engines you will find dozen of answers describing its use: https://ggplot2.tidyverse.org/reference/facet_grid.html and https://ggplot2.tidyverse.org/reference/facet_wrap.html
To prepare the faceting in ggplot2, a possible way is to bind both dataframes together by first specifying which group they belong such as:
cytokine$Group = "Cytokine"
il17$Group = "IL17"
DF <- rbind(cytokine,il17)
DF

    gene       variable value    Group
1:  CCT5     Res_immune 10.60 Cytokine
2:  CCT5 Non_res_immune 10.50 Cytokine
3:  CDK1     Res_immune  9.25 Cytokine
4: CXCL5 Non_res_immune  9.82 Cytokine
5:  GBP1     Res_immune  9.04 Cytokine
6:  GBP1 Non_res_immune  9.43 Cytokine
7:  CDK6     Res_immune 17.60     IL17
8:   IL6 Non_res_immune 14.50     IL17
9: CXCL8     Res_immune  1.25     IL17

Then, use facet_grid function of ggplot2:
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_grid(Group~., scales = "free")

Alternatively, you can use facet_grid (or facet_wrap) to have your plot in a single row:
library(ggplot2)
library(ggrepel)
ggplot(DF, aes(x = variable, y = value, fill = variable, label = gene))+
  geom_boxplot(alpha = 0.2)+
  geom_point(position = position_jitter(width = 0.2, seed = 2), show.legend = FALSE)+
  geom_text_repel(position = position_jitter(0.2, seed = 2), show.legend = FALSE)+
  facet_wrap(.~Group, scales = "free")

Reproducible data example
structure(list(gene = c("CCT5", "CCT5", "CDK1", "CXCL5", "GBP1",
"GBP1"), variable = c("Res_immune", "Non_res_immune", "Res_immune",
"Non_res_immune", "Res_immune", "Non_res_immune"), value = c(10.6,
10.5, 9.25, 9.82, 9.04, 9.43), Group = c("Cytokine", "Cytokine",
"Cytokine", "Cytokine", "Cytokine", "Cytokine")), row.names = c(NA,
-6L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

and
structure(list(gene = c("CDK6", "IL6", "CXCL8"), variable = c("Res_immune",
"Non_res_immune", "Res_immune"), value = c(17.6, 14.5, 1.25),
    Group = c("IL17", "IL17", "IL17")), row.names = c(NA, -3L
), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x55bfc0db9350>)

Sorry how I can place responders and non responders for each group next to each other? I mean for example for OBP responders and non responders being placed next to each other


text:
VCF's unknown ('.') values causing problems for GATK's Liftover Vcf
I am having a tough time calling LiftoverVcf in GATK. I have concatenated all of the Vcfs of interest and removed rows with problematic values but now when I run LiftoverVCF I get the following error:
Intel; Inflater: Intel; Provider GCS is available; Picard version: Version:4.0.12.0
INFO 2020-02-10 14:36:03 LiftoverVcf Loading up the target reference genome.
INFO 2020-02-10 14:36:16 LiftoverVcf Lifting variants over and sorting (not yet writing the output file.)

[Mon Feb 10 14:59:31 EST 2020] picard.vcf.LiftoverVcf done. Elapsed time: 23.50 minutes.
Runtime.totalMemory()=7414480896
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
java.lang.IllegalStateException: Key . found in VariantContext field INFO at 1:121387974 but this key isn't defined in the VCFHeader. We require all VCFs to have complete VCF headers by default.
at htsjdk.variant.vcf.VCFEncoder.fieldIsMissingFromHeaderError(VCFEncoder.java:202)
at htsjdk.variant.vcf.VCFEncoder.write(VCFEncoder.java:141)
at htsjdk.variant.variantcontext.writer.VCFWriter.add(VCFWriter.java:248)
at picard.vcf.LiftoverVcf.rejectVariant(LiftoverVcf.java:456)
at picard.vcf.LiftoverVcf.doWork(LiftoverVcf.java:357)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:295)
at org.broadinstitute.hellbender.cmdline.PicardCommandLineProgramExecutor.instanceMain(PicardCommandLineProgramExecutor.java:25)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

Not sure why I am getting this error other than there may be a . in one of the info fields. Here is what the first several rows of my VCF looks like minus the header fields:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT AltaiNea
1 10001 . T . 51.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=6.32;MQ0=242 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:21.05:0,21,265:1,2:0,1:1,0:116,129:2
1 10002 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.01:0,6,69:171,184:0,9:0,0:0,2:0
1 10003 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.02:0,6,72:222,201:1,0:0,0:0,8:0
1 10004 . C . 39.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.61;MQ0=247 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:9.02:0,9,97:0,1:268,211:0,0:0,1:0
1 10005 . C . 41.99 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.31;MQ0=246 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:12:0,12,119:0,0:341,287:0,0:0,0:0
1 10006 . C . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:6.02:0,6,72:0,1:408,347:0,0:2,0:0
1 10007 . T . 39.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.21;MQ0=246 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:9.02:0,9,97:0,0:0,0:0,0:506,451:0
1 10008 . A . 33.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.05;MQ0=247 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:3.01:0,3,40:547,507:0,0:0,0:0,0:0
1 10009 . A . 36.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.31;MQ0=248 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:6.02:0,6,70:600,534:0,1:0,0:0,1:0
1 10010 . C . 42.01 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=4.72;MQ0=244 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:12.02:0,12,125:0,0:649,547:0,1:1,0:0
1 10011 . C . 51.02 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=5.12;MQ0=242 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:21.03:0,21,209:0,2:732,612:0,1:0,0:0
1 10012 . C . 44.99 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=3.94;MQ0=243 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:14.99:0,15,141:0,3:804,657:0,0:6,1:0
1 10013 . T . 57.04 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.20;MQ0=238 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:27.04:0,27,285:0,0:1,0:0,0:904,760:0
1 10015 . A . 60.04 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.19;MQ0=238 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:30.04:0,30,301:985,850:0,0:1,0:0,0:0
1 10016 . C . 72.03 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.99;MQ0=235 GT:DP:GQ:PL:A:C:G:T:IR 0/0:249:42.04:0,42,427:2,1:1024,852:0,0:1,0:0
1 10017 . C . 84.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.37;MQ0=232 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:54.06:0,54,534:0,2:1100,912:0,0:0,0:0
1 10018 . C . 108.02 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.87;MQ0=224 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:78.03:0,78,757:0,0:1167,978:0,1:2,0:0
1 10019 . T . 96.05 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=10.24;MQ0=219 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:66.05:0,66,646:1,0:1,0:1,0:1258,1090:0
1 10020 . A . 65.97 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=9.74;MQ0=221 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:35.97:0,36,331:1309,1150:0,0:1,1:0,0:0
1 10021 . A . 65.34 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=8.14;MQ0=228 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:35.35:0,35,489:1339,1179:3,0:0,0:1,0:0
1 10022 . C . 92.98 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.79;MQ0=229 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:62.99:0,63,590:0,1:1389,1188:0,1:0,0:0
1 10023 . C . 98.93 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.69;MQ0=227 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:68.94:0,69,636:2,2:1453,1242:0,0:1,2:0
1 10024 . C . 98.95 . AC=0;AF=0.00;AN=2;DP=250;DS;MQ=7.46;MQ0=227 GT:DP:GQ:PL:A:C:G:T:IR 0/0:250:68.96:0,69,620:1,3:1497,1302:1,0:3,2:0

Is there a way for me to add . into the VCF header as a tag for all of the appropriate fields so that I can bypass this problem? I know from reading online that the . means that value is unknown. I would like to figure out how to instruct GATK to simply accept it regardless of the fact that it cannot be found in the VCF header.
Also, I am open to other ways I can liftover this genome to hg38. Below is the command I am using to do so using GATK.
java -jar ~/work/progs/gatk-4.0.12.0/gatk-package-4.0.12.0-local.jar LiftoverVcf -I AltaiNea.hg19_1000g.concat.LowQualRemoved.vcf -O AltaiNea.hg38_1000g.concat.LowQualRemoved.vcf -C hg19ToHg38.over.chain -R hg38.fa --REJECT rejectedRecords.vcf

You can easily get the header of the VCF printed into a text file using bcftools view -h file.vcf > new_header.txt , then manually add the INFO filed to the header following the required format (have a look at the VCF format specification), and then use bcftools reheader to get it replaced into the VCF file.
Providing the line at position 1:121387974 would be helpfull if you require help on how to build the new INFO line that you need in the new header

My problem was solved using --ALLOW-MISSING-FIELDS-IN-HEADER.


text:
Question: How to simulate 100k samples having 40 million SNPs in a proportion of case:control=30:70?
Note: this question can also be found on Biostars
I need to perform a stress test in a GWAS tool and the duty demands a dataset (plink format) having 100 thousand samples, having 40 million SNPs in a proportion of case:control=30:70. I'm performing the command:

plink1.9 --simulate ds1.sim --make-bed --out ds1 --simulate-ncases 30000 --simulate-ncontrols 70000

But after about 3 hours execution, I'm facing Error: File write failure., even though my environment has 256T of disk. (It's a container with a S3 (Object store) mounted through s3fs).
I was wondering if there's a clever way to do that? Perhaps breaking it into small chunks and merging up later somehow? It would prevent me to be waiting hours till an error happens and I loose either time and the sample generated so far.
Thanks in advance!

First comment: you should be using real data, rather than simulated data, as we are not very good at modelling genetic variation and recombination at a genome-wide scale. A way that this can be done with GCTA is mentioned in the PLINK manual.
I'll do some calculations to estimate the size of data you need to store:
Assume that this plink-formatted file is in tped format (one row per SNP), and the entire genome is diploid. That means that a single SNP/individual pair will take up three bytes. You want 30,000 cases and 70,000 controls, so each line from each sample will take up 300kB (plus a bit of change for the header). With 1 thousand SNPs, that's 300MB per sample; with 1 million SNPs, that's 300GB per sample; with 10 million SNPs, that's 3TB per sample; with 40 million SNPs, that's 12 TB for one sample.
Reading the plink manual, I see that the BED format is a two-bit compression for each SNP/individual pair, so one byte per 4 individuals, or about 12x compression (i.e. 1TB for one sample).
Okay, that's a lot. It's well within a 256TB storage limit, but might hit file size limits depending on the file system.
If there is a data generation issue with PLINK, it might be worth trying (as you suggest) splitting the input for simulation then merging afterwards. My guess is that doing the simulation on fewer SNPs (rather than fewer individuals) would work better. Plink has a --bmerge option for merging the compressed binary file formats.
I usually approach problems like these by working my way up from something that I know works. Try first merging two datasets of 1000 SNPs. If that works, try larger datasets. Rinse and repeat until the computer crashes, you get bored, or you reach your goal.

Hi @gringer thanks for the answer! Do you know where can I find free and open license real data with this amount of samples? I have been using the 1000 genome phase 3, but they have only 2504 samples. For this one, I managed to create synthetic phenotypes through some scripts I put together (https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb), but  GCTA --simu-cc/--simu-qt probably would address it in a more reliable way. Thanks for that!

For now, I managed to create a 20k samples with 8mi SNPs with Plink (1/5 from my final goal). My limitations is actually the S3 file size limitations. Also, your maths make sense, and suggests my scenarios doesn't. Thanks for that!
Btw - loved your TEDxWellington talk :)


text:
Algorithmic Information Measurement of Human Genome Databases?
What is the highest data compression ratio achieved on a population of genomes measured to include the compressing algorithm itself?  The size of the compressed database approximates its algorithmic information content*.
I've seen a claim of on the order of 1000:1 compression but it appears dependent on an external "reference" database.  In my opinion, a reference data must be included in the size of the compressing algorithm.
*The Universal Turing Machine notion of Minimum Description Length is also called the Kolmogorov Complexity or Algorithmic Information measure of the data.  I'm using the term "Algorithmic Information" as the most concise and accessible descriptive phrase for this measure.
The relevance of this measure is that according to the closely related theory of algorithmic probability, the algorithm that maximally compresses a dataset also provides the most optimal data-driven predictions.


text:
How to find a homolog that has a PDB structure available
I tried to find a homolog that has a PDB structure available so I can use this PDB file for comparative modeling. I ran a BLAST search but none of the search results seemed to have a known structure. Is there any way to find the results I want?

If you are using NCBI Blast and blastp against the pdb database gives you nothing, you can do a psiblast against the normal refseq database, do a few rounds, save the PSSM matrix and submit a new job with the PSSM matrix but against the PDB database.
You might need to limit the species of you get lots of too close homologues —say blasting E. coli protein for which you'd want to exclude Enterobacteriales.

You could use HMM-based comparison, for example with HHpred. This is generally the best way to find structural templates because structure is more conserved than sequence, and comparing HMMs will turn up hits of the same structure but minimal sequence overlap.
If HHpred turns up nothing you are in the region of de novo modelling, which is a different problem. However, it has advanced in recent years to the point where it may be useful to you. For example, you could try our method DMPfold.

To expand the first answer by @MatteoFurla: re-run the P-blast and under database option, which is set by default to non-redundant database (nr) select Protein databank (PDB), the third option from the bottom of the menu. All the Blast results will be from PDB and thus all results will have associated structures. Also note that the Genbank codes will be PDB codes.
If you have zero hits then consider PSI-Blast, how I suspect you are running P-Blast via the nr database.


text:
Vegetable databases usage
Good morning,
I'm aligning with blast some species, against all nt database to analyse the vegetables. However, in some cases we have species of animals mixed with vegetables and want only to detect the vegetable species. Normally, we would only expect to detect the vegetables only. But this is normal since we are using the all nt. So now, could anyone point me out, the way to detect only vegetables? Thanks for the help.

Cross posted at BioStars: https://www.biostars.org/p/421717/

vegetables are plants, meat are animals, you can filter nt to the taxonomic groups you are interested in

why was this question avaluated with -1?

-1 means it got downvoted, probably because a user did not appreciate the question, maybe because of lack of details, language, or simply because they found it is fun. I personally do not like the principle of downvoting questions at all, as it is utterly disrespectful. Still that is how SE works.
Speaking of your question, please consider being more precise in the future. Use correct terminology (plant, not vegetable, animal, not meat) and also say that you want to do this remotely, not via the web browser as the user who gave an answer assumed.

Hello Atpoint, thank you for the suggestions and I will do it.

Using the blast command line tool, you can limit your search to a specific taxon using the –taxids option.
So, in your case -taxids 3193 should limit your search to plants.

If you are using the ncbi blast webinterface, you can restrict your search by entering "plants" into the Organism field in the "Choose Search Set" section

Yes, but I want to use that restriction remotely...

Ah okay, I edited my answer to include the option in the command line tool

Thank you Wayne_Yus for the help! It was really helpful and I'm sorry if my question description mislead you.