text:
Feature extraction methods that can handle inconsistent numbers of atoms for molecular dynamics
I want to compare the protein dynamics
) pH 7 versus pH 3, or
) wild type versus mutant
The protein will have slightly different number of atoms at each condition, due to protonation or mutation, which is bad for machine learning. The alternative is to only consider the alpha-carbon atoms, but this will lose significant amount of information (e.g. side chain, salt bridge).
Is there a way in molecular dynamics to extract features from ensembles of different number of atoms?
Note: I have used common methods to decipher the difference between two conditions, including RMSD, RMSF, Rg, secondary structure, salt bridge, SASA, etc. They can explain the difference to some extent. However, I feel the selection of properties is "subjective", and there must be some hidden features not known. So I am seeking a machine learning method, that can "automatically" capture the most essential properties based on the coordinates as inputs.

This answer is no longer on topic. The question is about what is the best machine learning algorithm to use to analyse atomic coordinates from an MD trajectory in order to infer a novel properties.

Also, it should be preface for any confused reader that unfortunately "feature" is a technical term in both genetics and statistics. For the former a feature is a region of DNA, RNA or protein that has a given property, such as open reading frame or PFAM domain. While in the latter a feature is a preferably independent property that forms part of a feature vector for different types of machine learning.

Heavy atoms are common
Assuming the pH difference only results in protonations, the protons differ between the two conditions, but not the heavy atoms. These will retain all the side chain information.
Merit of backbone only
(Section no longer relevant due to OP edit)
The common metric used to compare two structures is root mean square deviation (RMSD) of the alpha carbons to assess how much does the backbone differ and with analysis of MD simulations it is used to assess whether one is more dynamic than the other (e.g. this Bioinformatics SE question). The reason why the RMSD is done with the alpha carbons and not the minimum common subset of atoms is that backbone changes are large scale, while the side chain changes only affect locally —so their inclusion would needlessly dilute the score. The C prime and backbone nitrogens behave like the C alpha so there is no need for these.
Extraction of coordinates
Assuming by features you mean atomic coordinates, there are loads of ways. Using Pymol module for Python, one could do something like:
import pymol2
with pymol2.PyMOL() as pymol: #can be parallelised
    pymol.cmd.load('trajectory.gro', 'mytraj')
    pymol.cmd.remove('h.') #remove hydrogens
    pymol.cmd.intra_fit("name CA", 1) #align all steps in the trajectory
    for i in range(1, pymol.cmd.count_states('mytraj')):
        for atom in pymol.cmd.get_model(f'state {i}').atom:
            atom.coord #List[int int int]

One could change pymol.cmd.get_model(f'state {i}') to select state {i} and name CA for Cαs only.
Edit: Mutated residues
In the case of mutated residues or residue differences due to the two proteins being close homologues, "mutating" the differing residues to alanine, i.e. by stripping all heavy atom beyond the Cβ (keeping C′, Cα, H, N, O & Cβ) in the MD trajectory —or glycine if one is glycine. If they are homologues with different lengths, this approach would mean skipping residues that match up with gaps in the MSA. Even though not uncommon for statistically analyses, it does pose some issues statistically. It may unbury hidden residues, but length variant generally happens with loops.
In the snippet above, the pymol mutagenesis wizard could be used, but that is overkill, does not work in parallel and the approach does not port to other PDB reading modules in whatever language one wants. So say you have a list to_be_trimmed of integers of PDB residue indices (not pose residue indices):
if atom.resi not in to_be_trimmed or atom.name.strip() in ('C', 'CA', 'H', 'N', 'O', 'CB'):
      do whatever with atom.coord

In acidic environments, you get some funky things happening like exposed aspartate-glycine residue pairs forming isoaspartate or exposed aspartate-proline resulting in a cleavage, etc. For these cases, which are nearly never modelled, ignoring these values is the best option.
Edit: small comment on SASA
The investigated property list is very extensive and some reader may adopt them. So I thought I should mention that per residue SASA is generally converted to relative surface area, which is residue independent, by dividing it by the maximum SASA for that residue (values from Tien et al. 2013: 'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}), where RSA > 0.2 is an exposed residue.

thanks, I have updated my question.


text:
Protein PDB/CIF text file reader that can fold a same section into one line
Is there a software that can recognise the protein PDB/CIF files, and fold multiple lines of a same section into one line, so that it is easier to handle? I need to understand the format in CIF so as to create my own CIF file.
For example, when Notepad++ is used to read a .py file, it can fold the same section into one line.
So for a PDB/CIF file, e.g. all the atom lines for a same chain can be folded into one line, and I can expand it if needed.
Like this picture, assuming "This PC" is the whole protein, while the other objects are the different chains.

You could try opening the PDBXML version of the file in an XML viewer that lets you collapse sections.

Thank you, the PDBXML file is 10 times the size of a CIF file. The 5TX1 has 100 MB. Can you recommend a software that is good for PDBXML?

Programs that work with PDB and mmCIF files may present the structure as a tree, for example:

(this is a Find dialog from ccp4MG)
but I'd be surprised if a text editor had such a feature. Normally, you don't edit PDB files in a text editor.

Thank you, I need to understand the format in CIF so as to create my own CIF file

@lanselibai do you want to create an mmCIF file in a text editor? why?

I am trying to do homology modelling for a virus capsid based on an existing CIF (5TX1, as the template) in the PDB database. The CIF template contains the assembly information, which I want to understand. But the CIF file is too long, which is very easy to get lost.

@lanselibai to understand assemblies you may start with Introduction to Biological Assemblies and the PDB Archive at RCSB's PDB101.

Sure, I will read it! Thank you!

source: Stack Exchange

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

sort: votes


text:
Custom-made kraken2 catalogue
I am looking for instructions on creating a custom-made kraken2 catalogue from the sequences of my choice and the corresponding phylogenetic tree in newick notation.
I will appreciate an input on whether/how it can be done. The prescription provided in the kraken2 manual seems to assume that sequences and phylogeny originate from the NCBI. Struo2 extends this to GTDB sequences. However, in my case the sequences are collections of concatenated genes (i.e., not full genomes), provided by a third party.
Remark: (in response to comments) I don't need to work with newick format - I have the tree in newick format. So converting it to dmp might be a part of the answer.

Newick has a lot of different ways to represent data, so it's helpful to know specifically what your input data looks like. When asking a question, it's helpful to provide both example input and example output. You've provided good information on the output side through the link, but no information on the input format.

@gringer Does kraken2 need this information? The question is not about solving my particular research problem (which is off topic?), but about creating a kraken2 catalogue from sequences and phylogeny, other than those readily provided by NCBI or GTDB (via Struo). [contd.]

You seem to be suggesting, that, once phylogeny is converted to dmp format, the problem is easily resolved - if that's the case, then state this explicitly, specify the tools for converting newick to dmp, and the minimal information that a newick tree should contain for creating a meaningful kraken2 catalogue. The would make for a nice answer, and useful for other users that stumble on the same problem.

I'm not trying to give answers, just clarify the question. My comments are more generally applicable to all questions on this site: a good question for StackExchange has a specific, well-defined input, and a well-defined output. All I'm trying to do is make sure that this question will receive good answers.

Yes, phylogeny is important. Kraken2 is a phylogenetic mapper, and uses the tree structure to work out the best match when it cannot be assigned reliably to a single species.

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


text:
Analyzing and visualizing metagenomic data
I am starting out in metagenomics,
Aim I wish to understand the whole analysis from the raw reads to interpretation.
Background My current pipeline is used for:

analyzing human gut microbiota cleans reads,
filters out human genome, and
maps the reads to known microbial species using IGC and mOTUs.

Question I am looking for basic tools to visualize the output, in order to get some biological and statistical insights (Linux and python are preferable environments).
For example, an obvious step would be to display the mOTUs analysis results on the phylogenetic tree, using the full phylogeny provided with this package. Whilst I could implement it myself, I was wondering if anyone else had previous experience?
Update
One package that I found useful is MOMR (through in R).


text:
Server for finding kmers in set of sequences
Is there a server/website somewhere where I can submit a list of DNA/RNA sequences and find the list of kmers hits and organisms where it's found? I checked the kraken website, but they don't have a webserver for it.
I am specifically looking for kmers, as a complement of just doing BLAST on the NCBI nt database or similar.
EDIT: So far I have seen:

the Kraken app on Illumina BaseSpace, but it's not trivial to upload FASTQ data that hasn't been produced from an Illumina instrument.
the CLARK trio of apps in InsideDNA, tried once with default parameters but produced empty results.

Why you need a webservice? Would not be easier to calculate them by your own?

One Codex is powerful and user friendly, and works in a similar way to Kraken underneath
https://www.onecodex.com/
The Kaiju web server might be of interest and is more sensitive than exclusively k-mer based methods
http://kaiju.binf.ku.dk/server
The title of this question is misleading – it sounds like you're actually looking for a webserver for taxonomic classification using k-mers.

Ah yes! I now remember having used onecodex in the past. I'll try both.

I tried onecodex both uploading a file and also copy+pasting the sequence into the 'Search' menu, and both came with no matches or 'Error'. Now trying kaiju...

I tried kaiju with standard parameters, no hits. I tried with relaxed parameters and got some hits. It did what it needed to do, and quite quickly as well. Thanks!

Glad that worked for you. k-mer based approaches have limited sensitivity in the presence of genetic variation and/or sequencing errors. Kaiju isn't exclusively k-mer based unlike One Codex IIRC. You could always try a more sensitive approach such as BLAST.

Perhaps you could edit the title to make this question more useful to others.

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


text:
What's a good ontology for drug names?
I have ...

A database with patient phenotypes in, stored as HPO terms
Genetic data in whatever format I need

I want ...

To store drug names in a way that won't make my life difficult
If it's feasible, to make them so that you could draw links between HPO terms, medications, genetic changes, and outcomes.  Note that I'm not actively trying to do that, but I don't want to exclude the possibility of it in the future

Ideally, I'd like those drug names in some kind of ontology, and even more ideally an ontology that's available on http://bioportal.bioontology.org
So - what's the best drug ontology to use in this case?
Update

What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease?
  Drugs by disease would be the most useful.  The information that would be most useful would be something like

Atenolol

Atenolol 25mg
Atenolol 50mg
Atenolol 100mg

Propanolol

Propanolol 10mg
Propanolol 40mg
Propanolol 80mg
...

In short, being able to get a list of dosages for a given drug.  It'd be extra-useful if I could also retrieve a set of drugs for a given condition, so in the above example to retrieve a list of beta blockers (Atenolol, Metoprolol, Nebivolol, ...)

Have you looked at any of the various drug ontologies listed in bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

MDDB doesn't seem to be structured as above
Likewise The Drug Ontology
RxNorm seems to be the best available that I've seen

What relationships would this drug ontology describe? Structural similarity of the compounds? Affected pathways? How would they be linked? Do you want all drugs that can cause the same side effect linked, for instance? Or all drugs that can be used for the same disease? Have you looked at any of the various drug ontologies listed in http://bioportal.bioontology.org? Which ones have you discarded so we don't suggest the same ones? What were they missing that you need?

@terdon - added some more information to the question, hopefully that answers your questions.  Apologies for the slight vagueness in some of my answers, the reason is purely because I don't know - for me it's almost just a popularity contest as to which is the best for me!

The thing is, we don't really do such broad questions on the Stack Exchange network. There's even a specific close reason for them ("too broad") so this is not the right place to ask for a general discussion of the pros and cons. We might be able to help you with a specific case, but not on "what's best" which is just a question of preference and will depend on what you want to do. What's best for you might not be best for me. So please make it as specific as you can.

@terdon - changed the wording, is that better? :)

Since this question's been asked there's been a (what looks like a first useful?) release of Drug Ontology DRON. Perhaps worth checking again?

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


text:
Adding matrix as a confounder in EnrichmentBrowser
Brain data, particualrly human brain data have very heterogeneus cell types, so it's important to normalise by cell type/add proportoins of different cell types to the formula when performing differential expression. We have some bulk RNA-seq brain data from disease samples and wildtypes. I have used the package BrainInABlender to get a matrix of the proportions of each cell type in each sample

Cell type
Sample_1
Sample 2

Neuron
0.5
0.7

Microglia
0.5
0.3

The rows are the types of cells, the columns are the samples. I wanted to know what's the best way to add this as a confounding variable in Enrichment Browser (or if that's not possible in DESeq2).
Edit: Current way I've tried to correct for this is in DESeq2. I have transposed the matrix of cell type proportoins then set DESeq's formula:
~ Condition + Condition:Neuron:Microglia:etc

This has drastically reduced the number of genes found to be differentially expressed but I don't know if it's the right way to do it.
Edit 2: When I tried removing the interaction terms, replacing them with +
~ Condition + Neuron + Microglia etc

As suggested on the Bioconductor support section the number of genes went up by an order of magnitude from just using condition which is the opposite of what should have happened. This persisted whether condition was at the start or end of my formula.

What have you tried? DESeq2 works with count data not decimal numbers. What is your biological question/goal?

I have no idea what is an enrichment browser. If you are referring to RNAseq analysis,  in DESeq2, when you set up the data with colData = ... just add the cell type neuron / microglia as a numeric column?

@llrs Sorry for not being clear eough. I have a seperate matrix of count data. I also have this matrix which has the proportoins of each cell type in each sample. My goal was to best correct for differences in proportoins of cell types in each sample.

@StupidWolf Enrichment Browser is an R package that allows for running several differential expression methods easily. I have added it to the coldata and seperated with : in my formula (Condition:Neuron:Microglia:etc but have no idea if that's right.

@Sethzard this doesn't solve the last question. Correct for differences in proportion of cell types for what? To see what cell type express?

@llrs Different types of cells express genes at different levels. Brain samples vary significantly in the proportions of cell types from one sample to the next. A sample containing a higher proportion of one cell type may report changes associated with differences in cell types rather than true differential expression unless samples are normalised or changes in those proportoins are corrected for. We have brain samples from a control group and a disease group and want to work out the truly differential expressed genes. I've edited the post to try to clarify too.


text:
CIBERSORT runtime error eval failed
I was running CIBERSORT but caught this error right after it started permutation:
org.rosuda.REngine.Rserve.RserveException: eval failed, request status: error code: 127
88

Ironically, I found that this problem was caused by the fact that the tables are too large!
They really should have designed the correct error messages!


text:
Can I get the graph generated by cellranger
I’ve run the cellranger analysis pipeline on single cell RNASeq datasets. I can import the matrix and graph-based clusters into R. Doing this I can optimise the dimension reduction and plot cells with coloured by clusters generated by cellranger. I wish to optimise the graph-based clustering as well.
Is it possible to either:

obtain the graph generated by cellranger from the files generated. (I can only find the final clusters in the analysis output)
reproduce generating a similar graph (used for clustering) from the output matrix (preferably with R or Python.

With either approach, I am looking for a list of edges for nearest neighbors in the gene-barcode UMI matrix to use as input for a clustering algorithm.

Could you provide some data as an example or pick the data from the vignette to answer you would work for you? What have you tried to get the graph or reproduce it?

I’ve tried several ways to reproduce this graph in R but I’m unsure that I understand well enough what cellranger is doing (based on Python code in the GitHub repo) to generate inputs for the Louvian algorithm. I’m now unsure that this is necessary since any cellranger run generates various output files for the analysis in a specific directory structure. I think this question requires someone with familiar with this software already and has run it themselves. I cannot provide the data that I’m working with currently and the files are too large even if I could.

Hi @TomKelly. I don't know if I understood exactly the question. However, you may find useful to use scanpy for importing the data, run the Louvain clustering, and plot the data with the clusters. That is a Python module.

@TomKelly Could you provide your attempts and why they fail? This would make it easier for me to avoid this pitfalls (if I try to answer your question)

I don't understand what plot you want. You know where the T-SNE coordinates are, right?

No, I can perform tSNE in R. I want the graph that the graph-based clusters are derived from. I want to annotate clusters into it. Can these be found in the output data or does it have to be computed from the matrix?

Cell Ranger
You can't download the tSNE coordinates for cells directly from the Analysis tab of the fancy, polished .html document that Cell Ranger produces. If you have access to the machine on which the pipeline was run, you can grab the intermediate files. Everything is well documented by 10X. They explain defaults etc. on that page. The tSNE coordinates are available at:
analysis/tsne/2_components/projection.csv

and the contents look like:
$ head -5 analysis/tsne/2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300

Seurat
And the second point: Cell Ranger is replicating the Seurat workflow. Seurat uses the Louvain method for community detection by default. Check out the intro tutorial and see if you can follow along. This tutorial works with 10X data, in fact the first function call in the tutorial reads in the 10X matrices. If you are already plotting in R you'll speed through this tutorial. Seurat is great to work with, you'll have tSNEs and clusters in no time. And there's so much more you can do from there!

Yes I’m familiar with the Seurat workflow. I don’t understand the Python source code well enough to know which parameters were used by cellranger.

If you're unwilling to use a different software package, why mention that as a possibility in your second point? The matrices are very lightly analyzed data. The unfiltered matrices are counts for all cell barcodes. The filtered matrices are count for cell arcodes determined to come from real cells. You still have to normalize and perform quality control filtering on this data efore it's ready for its first round of dimensional reduction and clustering.

If you want information about their algorithms, see here. There is also info on that page about the matrices. This should help with your understanding of how Cell Ranger is massaging the data.

I’m considering whether it is necessary to compute the graph again of it has already been done by cellranger. If it has not been saved by this program, it will be needed to compute it as you’ve suggested. In that case, more detailed information will be helpful. I have (and will continue) to do it this way but it would be reassuring to know I’m on the right track and not reinventing a graph structure that I already have saved to disk somewhere.

Ah, I was thinking you just just had the final product. If you have access to the machine where the pipeline was run you can get the tSNE coordinates.

Yes I ran the cellranger pipeline myself. However, I’m not looking for the tSNE coordinates. I’m looking for a list of edges used for the graph-based clustering (I can already find the clusters).

I've managed to get want I needed with the Seurat::FindClusters function from the object@snn output. I've detailed my process below.

Files in a cellranger run
Cellranger removes the intermediary files after the clusters have been computed. This can be observed in the following path from the cellranger project directory. This contains empty directories.
SC_RNA_AGGREGATOR_CS/SC_RNA_ANALYZER/RUN_GRAPH_CLUSTERING

However, the clusters generated by cellranger can be found in the file:
outs/analysis/clustering/graphclust/clusters.csv

Clusters can also be "exported" from outs/cloupe.cloupe file into a CSV using the Cloupe Browser. These CSVs can be imported and used for plotting in R or Python.
The graph structure used to compute these (with the Louvian algorithm) cannot be found in the files generated by cellranger.
Reproducing graph-based clustering
The Suerat R package performs graph-based clustering. The shared nearest neighbour (SNN) modularity optimisation uses the same parameters as the Python scripts called by cellranger (k=30). This generates an SNN matrix which can be used to generate a graph structure from an adjacency matrix. For a Seurat object test the SNN sparse matrix can be extract with:
FindClusters(test, dims.use = 1:20, algorithm = 1,
             resolution = 0.6, k.param = 30, save.SNN = TRUE,
             reduction.use = "pca", metric = "euclidean")
SNN_Matrix <- as.matrix(ceiling(as.matrix(test@snn)
graph_structure <- igraph::graph_from_adjacency(SNN_matrix)

FindClusters from the Seurat package requires a Seurat object that has already has a dimension reduction run on it such as RunPCA. The algorithm = 1 is the Louvian algorithm as performed by cellranger. The algorithm, reduction.use and metric can be changed as mentioned in the Suerat package documentation. Clustering algorithms not yet supported by Seurat can be run on the graph_structure resulting from extracting the SNN matrix. It is not necessary to compute the edges separately, the Seurat functions will compute and return these.
For example, the Leiden algorithm can be computed as follows:
devtools:install_github("TomKellyGenetics/leiden")
library(leiden)
clusters <- leiden(SNN_Matrix)


text:
Derive a GTF containing protein coding genes from a GTF file with Exons and CDS
Why I need a compatible file
I’m trying to run
velocyto with the R package to analyse RNA velocity (cell trajectories) with single cell RNASeq data. I have performed single cell analysis from 10x Genomics data using cellranger.
I have successfully aligned the reads to get loom files and imported these into R. I can get the velocity from these files by following the vignettes. However, I cannot reproduce the RNA velocity analysis based on “gene structure”.
I’m working with a different organism to the example (not human or mouse) so the annotation data provided in the example does not work. I have a GTF file for the latest annotation of this organism. However, it only contains “exon” and “CDS” as features. This appears to be the source of the problem. The “find.ip.sites” function requires a GTF with “features” = “gene” and one of the “attributes” to be “protein_coding”. These requirements are hard-coded into the velocyto.R function.
I have the following GTF files from the AtRTD2 dataset. The chromosome labels match the loom files in R.
Chr1    TAIR10  exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    Araport11   exon    6788    7069    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7157    7450    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7564    7649    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7762    7835    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7942    7987    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8236    8325    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8417    8464    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8571    8737    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7315    7450    .   -   1   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7564    7649    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7762    7835    .   -   2   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7942    7987    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 8236    8325    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";

The specifications for a GTF or GFF2 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

What I’m looking for
Is there a way to derive a compatible GTF file containing protein coding genes from the exons and CDS? I want to produce a GTF which contains genes in the features and protein_coding in the attributes. Is it possible to do this with existing tools or scripts?

What I’ve tried so far
I can modify the source code of the “find.ip.sites” function to run on my GTF file with these features missing. However, this requires running internal functions from the package written in Rcpp and means my workflow override future updates to the package. Running the rest of the vignette returns errors as no introns long enough have been identified (setting the thresholds lower is also incompatible with the general linear models called). Therefore I think it is better to generate a compatible GTF or GFF3 file rather than alerting the source code of the functions.
While intended for GTF files, the package functions can import GFF3 files, despite being described for GTF in the documentation. I’ve tried generating a GFF3 file using gffread from Cufflinks and replacing “feature” = “mRNA” with “gene” and adding “protein_coding”. This also returns errors when running the velocity algorithm. It does not work for either the GTF file used as an input for cellranger or the file generated by it. Is there a way to annotate protein coding genes and intron/exon boundaries based on a GTF file containing only exon and CDS annotations?
With cufflinks 2.2.1, a GFF3 was generated
gffread - E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gff
sed -i '/mRNA/s/gene_name=AT/gene_type=protein_coding;gene_name=AT/g' AtRTD2_19April2016.gff'
sed -i 's/mRNA/gene/g' AtRTD2_19April2016.gff

This is the GFF3 files that I've tried:
# gffread - E AtRTD2_19April2016.gtf -o-
##gff-version3
Chr1    TAIR10  gene    3631    5899    .   +   .   ID=AT1G01010.1;geneID=AT1G01010;gene_type=protein_coding;gene_name=AT1G01010
Chr1    TAIR10  exon    3631    3913    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    3996    4276    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4486    4605    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4706    5095    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5174    5326    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5439    5899    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3760    3913    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3996    4276    .   +   2   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4486    4605    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4706    5095    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5174    5326    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5439    5630    .   +   0   Parent=AT1G01010.1
Chr1    Araport11   gene    6788    8737    .   +   .   ID=AT1G01020_P2;geneID=AT1G01020;gene_type=protein_coding;gene_name=AT1G01020
Chr1    Araport11   exon    6788    7069    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7157    7450    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7564    7649    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7762    7835    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7942    7987    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8236    8325    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8417    8464    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8571    8737    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7315    7450    .   -   1   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7564    7649    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7762    7835    .   -   2   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7942    7987    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 8236    8325    .   -   0   Parent=AT1G01020_P2

The specifications for a GTF or GFF2 are:
The specifications for a GF3 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqid - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seq ID must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
type - type of feature. Must be a term or accession from the SOFA sequence ontology
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
phase - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attributes - A semicolon-separated list of tag-value pairs, providing additional information about each feature. Some of these tags are predefined, e.g. ID, Name, Alias, Parent - see the GFF documentation for more details.

It would probably help of you could [edit] your question and show us a few lines of each of the two GTF files so we can see what you mean more clearly.

You probably can use AWK for this, is the "gene" feature similar as "CDS" or is it from "first exon start" to "last exon end".

I think the question and the "What I'm looking for" are a little confusing. Specifically, it sounds like you want to generate a GTF in "What I'm looking for" but you want to get the amino acid sequence of proteins from a GTF in the title. Is there any way you can clarify or make the title and "What I'm looking for" consistent?

We have a complete GTF for Arabidopsis in Ensembl Plants.

Thanks @Emily_Ensembl, I was able to get velocyto to compute "global" RNA velocity based on this reference sequence. In case anyone else runs into this issue, the loom file needs to be re-computed with same GTF file with the shell utility to be compatible with the R commands.

In your case, I would definitely suggest following @Emily_Ensembl's advice and using the Arabidopsis GTF from Ensembl. But for future reference, if an Ensembl GTF wasn't available, you could build something like this using the gtf class from cgat
The cgat module and dependancies can be installed by following the instructions here. Specifically, you can download and run their installation script. This will set up a dedicated conda environment to run the versions of dependancies needed. Further documentation can be found here.
# download installation script:
curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/master/install.sh

# install the development version (recommended, no production version yet):
bash install.sh --devel

This requires python3 and anaconda to install. You may need to run this script in the conda environment:
# show available environments
conda info --envs
# activate by path
source activate /home/user/local/bin/conda-install/envs/cgat-a

Then run the follow as a python3 script. Copy the contents in a file such as convert_gtf.py. Then run it in the terminal on your gtf file:
python convert_gtf.py genes.gtf > genes_new.gtf

#python 3.6.4

import sys

from cgat import GTF
from cgatcore import iotools

for gene in GTF.flat_gene_iterator(
     GTF.iterator(iotools.open_file(sys.argv[1])), strict = False):
    gene_start = min (e.start for e in gene)
    gene_end = max(e.end for e in gene)

    if any(e.feature == "CDS" for e in gene):
        gene_type = "protein_coding"
    else:
        gene_type = "non_coding"

    for exon in gene:
        exon_line.setAttribute("gene_biotype", gene_type)

    gene_line = GTF.Entry().fromGTF(gene[0])
    gene_line.feature = "gene"
    gene_line.start = gene_start
    gene_line.end = gene_end

    # Its not clear what the transcript_id for a gene line is.
    # Technically speaking, Ensembl GTFs do not follow the GTF
    # standard as GTF entries must include transcript_ids, but
    # the ensembl ones don't.

    if hasattr(gene_line, "gene_name"):
        gene_line.attributes["gene_name"] = gene[0].gene_name
    else:
        gene_line.attributes["gene_name"] = gene[0].gene_id

    gene_line.attributes["gene_biotype"] = gene_type

    print (str(gene_line)+"\n")
    for exon in gene:
        print(str(exon_line)+"\n")

That’s great to know. Since this pertains to using tools in R, it would help to have instructions to install the library and run for those unfamiliar with Python.

I've tested this script and made minor changes. This is what worked for me (to covert the GTF: I'm not sure if it's compatible with velocyto yet).

Thanks for the fixes. The things have changed recently and I lose track of what version things are on. I believe the CGAT team is working on simplifying the installation. A couple of points: I'm not quite sure what you are trying to achieve with the if hasattr(gene_line, "gene_name) bit: presumably if it already has a gene_name attribute, you don't need to set one? Also note that your resulting lines arn't going to have a transcript_id. While this is also true of Ensembl GTFs, it does mean these GTFs don't mean the  GTF specification.

I had issues with missing gene names and often TAIR IDs are used instead of gene symbols. So of this might be specific to Arabidopsis. Most important is the fixes to the loop for exons as these lines weren’t defined.

Its not that I don't understand that you might want to add a gene_name to the attributes, its just that currently you test for the presence of a gene_name and then set it to its already currently defined value if it already has one. This seems redundant.


text:
Associating SNP and GENE
Assuming I have SNPs data using hg19, how can I know which SNP belongs into which Gene?
The data looks like:
chr10_103577643
chr10_124712463

and so on.
I want to add a column of Gene, which would tell to which Gene the SNP belongs.
The file is not a VCF file but a Txt file.

Please [edit] your question and show us an example of your data. What do you have? A list of rsIDs? A VCF file? Some other format? How many SNPs?

Just annotate the vcf file.

Use the Variant Effect Predictor(VEP)

Please don't just post links as answers. It would be great if you could add a short explanation of what VEP is, how to use it etc. Otherwise, this is really a comment and not an answer and would become useless in the (admittedly unlikely) case where the VEP were to be discontinued or its URL changed. For instance, pasting chr10_103577643 into VEP doesn't work, so you could explain what format the data would have to be in.

You can use snpEff or snpSIFT java annotator also for hg19.

Please don't just post links as answers. It would be great if you could add a short explanation of these tools do, how to use them etc. Otherwise, this is really a comment and not an answer and would become useless if the tools were to be discontinued or their URLs changed.


text:
Sampling haplotypes
I am trying to simulate different genome of peoples, I have data (VCF files) of various genes from the 1000K Gene project.
I want to simulate different whole genomes i.e generate a new population by combining real haplotypes I have. I am wondering what is the best way to approach the problem. That is what is an efficient method of generating realistic (Not just randomly select parts and recombining them) new genotypes based on the real genotypes I already have.
I am using Bioconductor packages VariantAnnotation to read and manipulate the VCF files and TxDb.Hsapiens.UCSC.hg19.knownGene to determine the Genes positions.
The data looks as the following:
> gene58@rowRanges
GRanges object with 91 ranges and 1 metadata column:
              seqnames    ranges strand | paramRangeID
                 <Rle> <IRanges>  <Rle> |     <factor>
  rs551585351        1 229566998      * |         <NA>
  rs528384854        1 229567027      * |         <NA>
  rs542093083        1 229567063      * |         <NA>
  rs561849701        1 229567128      * |         <NA>
  rs531042647        1 229567160      * |         <NA>
          ...      ...       ...    ... .          ...
  rs565479298        1 229569784      * |         <NA>
  rs572772527        1 229569785      * |         <NA>
     rs605430        1 229569803      * |         <NA>
     rs605428        1 229569804      * |         <NA>
  rs368699658        1 229569810      * |         <NA>
  -------
  seqinfo: 86 sequences from "hg19" genome

The genotype matrix of the gene:
> gene.mat[1:5, 1:5]
            HG00867 HG02371 HG00759 HG00766 HG00844
rs551585351 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs528384854 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs542093083 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs561849701 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"
rs531042647 "0|0"   "0|0"   "0|0"   "0|0"   "0|0"

For simple haplotype simulation of unrelated individuals, there is Montana's venerable HapSim. The model used in that tool is simple and may be good enough for your purposes. You will need to acquire a recombination rate map from for example the HapMap consortium or from 1000G. A cursory google search also revealed this tool and walkthrough that may be suited to your problem: https://adimitromanolakis.github.io/sim1000G/inst/doc/SimulatingFamilyData.html
There are lots of other tools with more advanced models, but it sounds like those would be overkill here.

The data in your matrix appears to be from a VCF file (For reference, see https://github.com/samtools/hts-specs/blob/master/VCFv4.3.pdf).
In the VCF specification, a heterozygous genotype can be specified as 0|1. The zero and one indicate the allele numbers (0=ref, 1=first alt) that make up the genotype. The "pipe" (|) is used to indicate that the alleles are in phase with the previous entry in the file. Therefore, if you observe data like this:
"0|1"    "0|0"
"0|1"    "1|1"
"0|1"    "1|1"

Then this means that the first sample has the haplotypes [0,0,0] and [1,1,1]; and that the second sample has the haplotypes [0,1,1] and [0,1,1].

Thanks, but its not what I meant, I wanted to know how to recombine the haplotypes in order to generate new population out of it. Is there a place where recombination is more likely occur? Is there a program where I can input real haplotypes and it will generate a population?

Oh! Well there certainly are programs for that. Any chance you could rephrase your post to contain that question?

Hopes its better, let me know if I left anything unclear.

Ensembl have haplotypes from 1000 Genomes available by transcript.


text:
size of the pathways for analysis and filtration
I have recently started working on a substance's effect on a cell line in different dosages. for this, there is a tool called bmdexpress2 that I am using. Its input is the normalized counts from RNASeq for each dosage as a big matrix. When it comes to the pathway analysis step, unlike hallmarks of GSEA, this tool uses some databases for defining pathways which involve pathways of even 2-3 genes.
So what I want to discuss here is the strength of small pathways and filtering thresholds. How can we bioinformatically decide on the importance of a pathway made of 3 genes? Should we just filter them out as there are much more comprehensive pathways? Or do their GO Levels matter? Also there are many cases that out of those 3 genes, 1 of them is differentially expressed; is having 1 gene diff. expressed out of X genes in a pathway enough to say that the pathway is enriched?
Here is an example GO term for a small pathway.
And here is an example of the output I get from the tool. all genes platform column is the number of genes involved in that pathway, input genes column is how many of the genes in that pathway are significantly differentiated. I believe the rest is self explanatory.

Looking at the table, can we really say, ESCRT III Complex, for example is a good pathway to focus while having 1 gene out of 13 in a pathway but that change is significant?
Thanks

Why is the GO term relevant? Is that one of the terms of one of the proteins in your pathway? Is it a term that was reported enriched for that pathway? What sort of pathway analysis is this? Does the tool report over/under represented BP GO terms? In what? Does it return a list of pathways with over/under represented diff expressed genes? Presumably, the statistical test will only show something if it is significant, which will be hard with very small pathways, but please [edit] your question and give us some more details about the analysis.

statistically. the tool only returns fisher exact two tail test value and I already filtered my list of pathways according to that. the tool also returns pathways as go terms and the example go term is one of the significantly enriched

There's nothing special about the term. In order to help you, we would need to see the relevant pathway and the numbers involved. As I said before, please [edit] your question and give us more details.

Are you talking about pathway or about genes related to a GO term? It is not the same as pathway, although if you give us your definition of pathway it might be that it fits... If you don't take into account the structure of GO terms you'll get many false positives. Also, could you [edit] your question and include the data in text format? It will be easier to search more information about the genes and GO terms.. Also did you correct by multiple comparisons (as an answer suggest)?

The problem here is that any one GO term with a small number of members is unlikely to include any differentially expressed genes by chance, but if you test a large number of GO terms, some of them definitley will. The solution is to correct the p-values for multiple comparisons. If you do this, then none of those GO terms with a single DE member will be significant. The larger GO terms with multiple DE members in the chart you show probably won't survive correction either. Hopefully you have some that are more significant.


text:
How to plot genomic.fna fasta file in R using Gviz?
I downloaded the genome of Staph aureus as DNA sequence from NCBI. I would like to visualize it using Gviz. Is Gviz the right tool to visualize genomes? What are the necessary steps from the sequence file to an object that can be used in Gviz? I searched the user guide and the reference manual. It is stated that it is possible to use fasta file, but not how to use them.
We will generate lots of genomes that we will have to annotate and compare. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genomes. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something and get going. It is my first steps with genomic data, I don't even know what makes sense.
Here is the head of the file:
>NC_007795.1 Staphylococcus aureus subsp. aureus NCTC 8325 chromosome, complete genome
CGATTAAAGATAGAAATACACGATGCGAGCAATCAAATTTCATAACATCACCATGAGTTTGGTCCGAAGCATGAGTGTTT
ACAATGTTCGAACACCTTATACAGTTCTTATACATACTTTATAAATTATTTCCCAAACTGTTTTGATACACTCACTAACA
GATACTCTATAGAAGGAAAAGTTATCCACTTATGCACATTTATAGTTTTCAGAATTGTGGATAATTAGAAATTACACACA
AAGTTATACTATTTTTAGCAACATATTCACAGGTATTTGACATATAGAGAACTGAAAAAGTATAATTGTGTGGATAAGTC
GTCCAACTCATGATTTTATAAGGATTTATTTATTGATTTTTACATAAAAATACTGTGCATAACTAATAAGCAAGATAAAG
TTATCCACCGATTGTTATTAACTTGTGGATAATTATTAACATGGTGTGTTTAGAAGTTATCCACGGCTGTTATTTTTGTG
TATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCA
AGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAAAGATACTGAGCTTTACACGATTAAAGATGGTGAAGCTATCG
TATTATCGAGTATTCCTTTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGATGTT

what do you want to visualize? the proteins, the GC content, some other features etc.?

Ultimately, I would like to annotate protein coding regions and compare them to other genomes.

What have you tried to do? Which code have you tried and what errors and message you got?

I would like to visualize annotated genes and compare a set of genomes.

This question seems too broad. Could you expand what do you want to compare in a set of genomes, or what type of annotation do you want to see (coding regions or proteins or ORF?)? Also please decide if your question is about seeing annotations OR compare a set of genomes (to compare genomes one would usually calculate some metrics in each genome and compare them). Remember you can [edit] the question.

It is not broad. It is very specific. I am looking for an example to plot SOMETHING with Gviz. An example, for someone who has never used it before.

I am not familiar with Gviz, but I am familiar with a variety of tools for visualizing genome features. At the scale of entire microbial genomes, the primary sequence itself is an overwhelming amount of information and can't really be visualized in a meaningful way by itself. I mean, you could color the nucleotides, but few patterns in DNA can be discerned by an unassisted human eye.
What we're typically interested in are features encoded in the genome sequence: protein-coding genes, non-coding RNAs, transposable elements, regulatory elements, binding sites, epigenetic modifications, and so on. A particular scientist is usually only interested in a subset of these features, based on what particular questions they are investigating. In your comment you say you "would like to annotate protein coding regions and compare them to other genomes", so it seems like we're on the same page here.
Visualization tools often display the genome linearly (even, for example, if it's a circular bacterial chromosome), zoomed out so that the nucleotide sequence isn't visible, and overlayed with boxes, arrows, or other glyphs representing the genes or other features of interest. Here are some examples:

This question thread contains images generated with the GenomeTools AnnotationSketch program (gt sketch command).
The UCSC genome browser which provides web-based access to dynamic displays of the human genome and many other mammalian genomes.
The IGV browser can be downloaded and easily run on your personal laptop or desktop computer.
I'm not familiar with Gviz, but it definitely looks like it belongs in this class of tools.

In each of these cases, the visualization software is not doing the work to determine the positions of the genes and other genome features. This was done independently using tools designed for this task, storing the information in an annotation file in GFF3, GTF, or BED format. The visualization tools need both the genome sequence (in Fasta format) and one or more annotation files (in GFF3, GTF, or BED format) to make meaningful visualizations.

I am touching on genomics in a project. So I am quite new to the field and don't really know what is possible, feasible. Thanks a lot for your answer. I am working with bacterial genomes.

I use prokka to annotate the genome.

Prokka should create a GFF3 file, which you can use with the Fasta file to visualize your genome. But of course, you're not the first person to annotate S. aureus, so you may also consider downloading a reference annotation for the species from NCBI RefSeq.

We will generate lots of genomes that we will have to annotate though. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genome. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something.

Great. That type of background is very useful to include in the original question, so that the community knows how to focus their answers. I hope this helped!

I added it to the main question. Thanks for your comment.


text:
How to export peptide-shaker data from commandline?
I would like to export the data from Protein-Shaker proteomics analysis software and use it in downstream analysis. For example the overview data, numbers of validated peptides per protein, confidence, etc. Ideally I would like to generate this output from commandline without the gui and, if possible to plain txt or csv.

You have to start PeptideShaker calling java with the -cp argument an input project file and provide a comma/space-separated list of integers for the reports you want to export:
java -cp PeptideShaker-1.16.37.jar eu.isas.peptideshaker.cmd.ReportCLI -in input.cpsx -reports "0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10"`

Here is the list of reports you can export:
-reports                   Comma separated list of types of report to export.
                           0: Certificate of Analysis,
                           1: Default Hierarchical Report,
                           2: Default PSM Phosphorylation Report,
                           3: Default PSM Report,
                           4: Default PSM Report with non-validated matches,
                           5: Default Peptide Phosphorylation Report,
                           6: Default Peptide Report,
                           7: Default Peptide Report with non-validated matches,
                           8: Default Protein Phosphorylation Report,
                           9: Default Protein Report,
                           10: Default Protein Report with non-validated matches,
                           11: Extended PSM Report,
                           12-n: Your own custom reports.

I don't know whether you can also use the SearchGUI output zip to generate this files. The program accepts zip files as input, but I think that means PeptideShaker project that are exported to zip files, but not SearchGUI output zips. So, in the end you might need to run the GUI at least one time, unfortunately.


text:
Search ChebiEntity search by common or synonym name?
I can download information from ChEBI if I know the corresponding ChEBI ID:
from libchebipy import ChebiEntity
chebi_entity = ChebiEntity("15903")

And then look at all the synonyms that the compound has:
chebi_entity.get_names()
> [{'_Name__name': 'beta-D-Glucose', '_Name__typ': 'SYNONYM', '_Name__source': 'KEGG COMPOUND', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucose', '_Name__typ': 'NAME', '_Name__source': 'UniProt', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'beta-D-glucopyranose', '_Name__typ': 'IUPAC NAME', '_Name__source': 'IUPAC', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'BETA-D-GLUCOSE', '_Name__typ': 'SYNONYM', '_Name__source': 'PDBeChem', '_Name__adapted': False, '_Name__language': 'en'},
 {'_Name__name': 'WURCS=2.0/1,1,0/[a2122h-1b_1-5]/1/', '_Name__typ': 'SYNONYM', '_Name__source': 'GlyTouCan', '_Name__adapted': False, '_Name__language': 'en'}]

Is there a way to get the ChEBI ID with one of  the synonyms?
beta-D-glucose -> 15903

Is there any other python tool or web-api that can achieve this?

How very interesting ...
Yes its easy just download the flat file database here:

https://ftp.ebi.ac.uk/pub/databases/chebi/Flat_file_tab_delimited/names.tsv.gz

WARNING: it's a 5.5M file ... don't click the link unless you really want to download a file this size.
You can simply load it straight into a pandas dataframe from read_csv() and then interrogation it via loc. You could place the header NAME as the index. There are other way in, you can shift the dataframe to a dictionary and interrogate that way, or place it straight it SQL via e.g. sqlite3
Here I just loaded it onto Excel and found the entry you were looking for ...

....

Edited: sorry a pandas iloc() solution wasn't appropriate.


text:
How to run proteowizzard on linux using wine?
I am trying to run proteowizzard using wine to be able to convert raw mass spectrometry files to myML format. Following, these instructions` which seem terribly outdated. I am getting many issues. I am running Ubuntu 16.04.
We want to convert output from a ThermoFisher instrument to mzML or mzXML format.
> ~/software/proteowizzard>wine msconvert.exe
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:process:SetProcessShutdownParameters (00000380, 00000000):
> partial stub. err:module:import_dll Library VCOMP110.DLL (which is
> needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\baf2sql_c.dll") not
> found err:module:import_dll Library baf2sql_c.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library VCRUNTIME140.dll (which is needed
> by L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-convert-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-environment-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-filesystem-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library VCRUNTIME140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-filesystem-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library
> api-ms-win-crt-convert-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-environment-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library msparser.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:LdrInitializeThunk Main exe initialization for
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe" failed,
> status c0000135

Which issues did you find? You mention wine, so I assume you are using a linux computer to run a program for windows. Could you add more information about the requirements of this program?

It might be better to ask this on [ubuntu.se], but either way please [edit] and show us exactly what commands you ran and the specific errors you are getting. We can't help you if you don't tell us what the symptoms are.

I am wondering whether there is anyone here who runs proteowizzard on linux, who knows how to set it up.

Try this Docker container.
All the vendor formats work except T2D, YEP, and FID.

I sympathize ...
If it is a Mac just run 'PlayonMac' its a brilliant frontend for Wine without all the bugs and hassles it causes. Its free and available here


text:
Calculate molecular mass from chemical formula C9H15N4O8P with Python?
Can I calculate the molecular mass from a formula like this C9H15N4O8P with Python e.g. using rdkit? I think rdkit.Chem.Descriptors.MolWt needs a molecule object. Though, to generate that I would need e.g. a SMILES string. Is there any other python module that can do the trick?
I found this: http://verahill.blogspot.com/2012/10/molecular-weight-calculator-in-python.html

I tried to adapt this approach within rdKit to calculate exact mass (for mass spectrometry), but was not successful.
This library (https://pypi.org/project/molmass/) appears to be able to do this:
from molmass import Formula
f = Formula('C7H11N3O2')
f.formula
f.isotope.mass
'169.08512660696'

be cautious with molmass because it admits that it does not include the mass defect in the calculations so for large molecules there will be a sizable inaccuracy.

As you mentioned, rdKit does not accept a chemical formula as input. The main reason for this are probably isomers. You could however use rdKit to initiate all the atoms, get their mass and add them together.
import re
from rdkit import Chem

def get_mass(formula):

    parts = re.findall("[A-Z][a-z]?|[0-9]+", formula)
    mass = 0

    for index in range(len(parts)):
        if parts[index].isnumeric():
            continue

        atom = Chem.Atom(parts[index])
        multiplier = int(parts[index + 1]) if len(parts) > index + 1 and parts[index + 1].isnumeric() else 1
        mass += atom.GetMass() * multiplier

    return mass

Then...
>>> get_mass("C9H15N4O8P")
338.21299999999997


text:
How to convert a SDF molecular structure file to a pandas dataframe with Python?
What is the best python tool to convert a SDF file to a structured data frame (pandas.DataFrame). Similar to Molconvert from ChemAxon? SDF can have different datafields for each compound which I would like to have in a dataframe. The 2D or 3D structures can be disregarded.

Please provide a small example file

It is easily done with RDKit:
from rdkit.Chem import PandasTools

fn = '/mypath/dataset.sdf'

df = PandasTools.LoadSDF(fn, embedProps=True, molColName=None, smilesName='smiles')


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

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

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

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

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


text:
How to retrieve dataset names from python statsmodels?
The Python statsmodels library grants access to many R-datasets which can be downloaded e.g. with:
import statsmodels.api as sm
df = sm.datasets.get_rdataset(dataname='iris', package='datasets').data

But for this you have to know the id or the name of the dataset.
sm.datasets.get_rdataset??
Signature: sm.datasets.get_rdataset(dataname, package='datasets', cache=False)
Source:
def get_rdataset(dataname, package="datasets", cache=False):
    """download and return R dataset

    Parameters
    ----------
    dataname : str
        The name of the dataset you want to download

There is a list of available datasets here, but this list does not contain the identifiers. Also not the detailed view.
Is there an overview of all the available datasets that contains the ids e.g. "iris"?
There is an index that is supposed to list all the datasets, but e.g. 'rice' as a name does not work.

Here is the website of the creator of rdatasets I guess with the info you are looking for:
http://vincentarelbundock.github.io/Rdatasets/datasets.html
Here is a short look at it:
Package     Item    Title   Rows    Cols    n_binary    n_character     n_factor    n_logical   n_numeric   CSV     Doc
boot    acme    Monthly Excess Returns  60  3   0   1   0   0   2   CSV     DOC
boot    aids    Delay in AIDS Reporting in England and Wales    570     6   1   0   0   0   6   CSV     DOC
boot    aircondit   Failures of Air-conditioning Equipment  12  1   0   0   0   0   1   CSV     DOC
boot    aircondit7  Failures of Air-conditioning Equipment  24  1   0   0   0   0   1   CSV     DOC
boot    amis    Car Speeding and Warning Signs  8437    4   1   0   0   0   4   CSV     DOC
boot    aml     Remission Times for Acute Myelogenous Leukaemia     23  3   2   0   0   0   3   CSV     DOC
...

I assume the Item column are the identifiers.

As stated above, that does not seem to work. Try for example 'rice'.

Did you change the package parameter to DAAG for the dataset rice? Since it’s not in the package datasets according to the index site.

You can view all the available datasets at this URL https://github.com/vincentarelbundock/Rdatasets/blob/master/datasets.csv
The available ones in statsmodels are those with the tag datasets which looks like line 372 - 455.
I didn't test them all, but everyone I tried from a small random sample seemed to work fine.

Please, try this code to get all modules in datasets:
import statsmodels.api as sm
import inspect

def get_list_of_modules_in_datasets():
  # get list of strings for the objects in datasets
  list_str_in_datasets = sm.datasets.__all__
  print(list_str_in_datasets)
  # get list of objects in datasets
  list_objs_in_datasets = [eval('sm.datasets.' + dsstr) for dsstr in list_str_in_datasets]
  print(list_objs_in_datasets)
  # get list of modules in datasets
  list_modules_in_datasets = []
  for dso in list_objs_in_datasets:
    if (inspect.ismodule(dso)):
      list_modules_in_datasets.append(dso)
  print(list_modules_in_datasets)
  print("Number of Modules in datasets = %d" %(len(list_modules_in_datasets)))
  return (list_modules_in_datasets)

if __name__ == '__main__':

  list_modules_in_datasets = get_list_of_modules_in_datasets()

  # print data for all dataframes from the datasets
  for ds in list_modules_in_datasets:
    df = ds.load_pandas().data
    print(list(df.columns))
    print(df.shape)
    print(df)


text:
What is the best way to drop alignment artifacts in sequence objects in biopython?
I have a bunch of sequences created with Roary which I want to preprocess. When I read the sequences (genomics data) I want to translate them to protein sequences using:
import pandas as pd
from Bio.Seq import Seq
from Bio import SeqIO

def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.translate())} for rec in fasta_sequences])
    return sequences

The problem is that .translate() fails with:

TranslationError: Codon 'T--' is invalid

Is there a build-in function to ignore the dashes?

That is some cool comprehension within a dataframe ... of course you need import pandas as pd

I think the .ungap() is what I was looking for.
def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.ungap().translate())} for rec in fasta_sequences])
    return sequences


text:
How to create Phylogenetic Trees from fasta files in Python or R?
I have around a hundred Fasta files (and will collect several thousand) with DNA sequences and +50x coverage. What is a recommended method to construct a phylogenetic tree? Solutions in Python or R are sought.
I found Phylo from Biopython only handles already calculated trees.

Would you be open to using a program designed to do this or do you need to roll your own in a programming language?

I would not look for a package for this, but instead build a small pipeline calling external tools with something like the following workflow:

Cluster the ~100 sequences with CD-HIT-EST/PSI-CD-HIT or many other options
Take all the sequences that form one individual cluster and build a multiple sequence alignment (MSA) with MAFFT/ClustalOmega or similar
Take the MSA and build a phylogenetic tree with a Maximum-Likelihood approach like IQ-TREE or similar
Visualize the tree file with Jalview or similar

Of course this is rather general and depending on exactly what you're doing you may want a different workflow and/or different tools. You should also explore the parameter space, do not assume the defaults are necessarily good choices

Yes, I agree with this and thats what I do. R is has its uses however in phylogeny

Jalview can also make the tree for you. Just load the alignment into jalview and tell it to build a tree from it.

@terdon not recommended though, I think jalview only implements NJ

Oh sure, I am not saying it's the best tool. I just know that ~10 years ago when I last used it, it was perfectly capable of producing a decent tree. But then I was using about a dozen sequences, all of which were homologous. So not a very hard tree to build.

@terdon 10 years ago is about right - pre-history and even Perl programmers existed then ... I think you critiqued me for 'samtools' so touche!

I agree with the first two steps. For the tree building I might consider RAxML. I might then load the finished tree into R with the ape/phangorn package for manipulation/visualisation.  For stand-along visualisation figtree is great.

I agree with Chris Rands that a reasonable approach would be to call external tools.
However, if you really want to do the phylogeny from within Python, you could use the P4 package, which is a bit complicated to handle but gives you lots of options in the way to build MCMC-based bayesian phylogenies:
https://github.com/pgfoster/p4-phylogenetics
You would still need something else to align the sequences before.
To visualize the tree using python, you could use the ete toolkit, which is likely more powerful than what you can find in Biopython: http://etetoolkit.org/

The obvious single answer is R "ape". This will give you access to PhylML for tree building and Clustal/Muscle for alignment building. The paths to the binarys are important. There are several distance methods in there such as NJ and BIONJ. Its distance approaches however don't look mainstream, but I could be wrong.
There are functions within ape which are cool, the tree sorting is very cool and I need to read through this with much greater care. Personally I wouldn't perform a core phylogenetic analysis within R, because the standalones are sufficient and the analysis is intensive.
https://cran.r-project.org/web/packages/ape/ape.pdf

(if I understand your situation correctly) https://www.rdocumentation.org/packages/seqinr/versions/3.6-1/topics/read.alignment shows how to use the function read.alignment which can take fasta msf etc. The docs provide the example' read.alignment(file = system.file("sequences/LTPs128_SSU_aligned_First_Two.fasta",  package = "seqinr"), format = "fasta", whole.header = TRUE) but you can use this code below (assumes those files are aligned) to go from reading the tree to getting the distances, producing the neighbor joining phylogenetic tree, and then plotting the tree.
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

fasta.res <- read.alignment(file = "geneticAlignment.msf", format = "fasta")
fasta.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
fasta.res.dist.alignment.nj = nj(fasta.res.dist.alignment)
plot(fasta.res.dist.alignment.nj, main = "from fasta files")

Within the R/Bioconductor ecosystem you could use the DECIPHER package (https://bioconductor.org/packages/release/bioc/html/DECIPHER.html) to align your sequences and then continue with either ape or phangorn  (https://klausvigo.github.io/phangorn/) packages to compute a tree. You can then visualize your tree using for example the ggtree package (https://bioconductor.org/packages/release/bioc/html/ggtree.html).


text:
How to do protein quantification with Peptideshaker output (TMT labeled, MS2)?
I have TMT tagged MS2 mass spec run processed with Peptideshaker. I saved the project in a peptideshaker format cpsx. Now, I would like to do the quantification with Reporter. I can load the file with Reporter GUI. Though Reporter is still in development and cannot be run from commandline yet. Are there other tools available that can run from Peptideshaker exported reports for downstream analysis? I tried to run Reporter from commandline version 0.7.16. I will try 0.7.17 which was just released.

It seems to me there should be something, but without knowing what you have tried we could repeat the same search you did. What search did you perform? Do you need it to be accessible via command line?


text:
Snakemake fails because it cannot create logs inside a read-only Singularity Container
I am using Snakemake from inside a Singularity container. Because of that, the directory that Snakemake wants to log to cannot be created, and I get an OSError [Errno 30] Read-only file system: '.snakemake, more specifically the traceback
Traceback (most recent call last):
  File "/usr/local/bin/snakemake", line 8, in <module>
    sys.exit(main())
  File "/usr/local/lib/python3.8/dist-packages/snakemake/__init__.py", line 2562, in main
    success = snakemake(
  File "/usr/local/lib/python3.8/dist-packages/snakemake/__init__.py", line 525, in snakemake
    logger.setup_logfile()
  File "/usr/local/lib/python3.8/dist-packages/snakemake/logging.py", line 137, in setup_logfile
    os.makedirs(os.path.join(".snakemake", "log"), exist_ok=True)
  File "/usr/lib/python3.8/os.py", line 213, in makedirs
    makedirs(head, exist_ok=exist_ok)
  File "/usr/lib/python3.8/os.py", line 223, in makedirs
    mkdir(name, mode)
OSError: [Errno 30] Read-only file system: '.snakemake'

The Snakefile I am using is inside a Singularity-container (there's no way around that, I am on a hpc and don't have root, also using a --sandboxed container doesn't work), and I call it like this (I've also played around with binding other directories without success):
singularity run --bind $(realpath .) --bind ~/container.sif ~/container.sif test_simple --cores 4

The test_simple rule is a really simple one that only unzips a gzip, for testing - the problem however lies generally in Snakemake wanting to create the log at <snakemake_base>/.snakemake/log (see https://github.com/snakemake/snakemake/blob/bd5a9239b2684b57ab0143bb709e77afa1493cdc/snakemake/logging.py#L308), which in my case is inside the container.
It doesn't help to use an overlay --overlay overlay.img for some reason, and I've also tried using my own (empty) log-handle (--log-handler-script ~/loghandler.py) without any success.
Any ideas on how to circumvent this or disable the logging or write it to somewhere else?

I don't know what happens in your case, but I use Singularity containers with snakemake inside, and when it runs, the .snakemake folder is created inside the current working directory, not in a read-only folder. It might help seeing the container definition file.

Basic sanity check: do you run the container from a folder in which you have write access?

Can you write to the logs into /tmp/ ?

Or maybe don't run Snakemake in the container (e.g., conda install -c conda-forge -c bioconda snakemake doesn't need root).


text:
Mapping ICD9 codes in UKBB data to PheCodes using PheWAS package
I am working on the UKBB data and trying to specify cases and controls by converting ICD9 codes to PheCodes. I have input data (icdt) where the first column (f.eid) is patient IDs and the remaining columns are ICD9 codes. I am using CreatePhewasTable function in R for the conversion. The code seems working well, But I am not getting enough 'TRUE' values for the PheCode 250.1 and 250.2, which represent PheCode for T1D and T2D. I am not sure what I am overlooking because I am getting most of the 'FALSE' values for T1D and T2D columns (250.1 and 250.2). I am sharing my R script
library(data.table)
library(dplyr)
library(PheWAS)
icdt <- fread("icd_codes.tab",header=T,sep="\t")
icdt <- as.data.frame(icdt)

for (i in 215:261) {
  icd9 <- data.frame(matrix(nrow = dim(icdt)[1], ncol = 2))
  icd9[,1] <- icdt[,1]
  colnames(icd9)[1] <- 'f.eid'
  colnames(icd9)[2] <- 'DiagnosisCodes'
  icd9[,2] <- icdt[,i]
  icd9 <- icd9[!rowSums((is.na(icd9))),]
  icd.codes <- icd9 %>% mutate(DiagnosisCodes = (paste0(substr(DiagnosisCodes,1,3),".",substr(DiagnosisCodes,4,nchar(DiagnosisCodes)))))
  icd.count <-  as.data.frame(count(icd.codes, f.eid, DiagnosisCodes))
  phenotypes <- createPhewasTable(icd.count, min.code.count = 1, add.exclusions = T, translate=T)
  phenotypes <- as.data.frame(phenotypes)
  assign(paste('Phecodes',i,sep=''),phenotypes)
  myfile <- paste0( "ALL_Phecodes", i, ".txt")
  write.table(phenotypes, file = myfile, sep = "\t", row.names = FALSE, col.names = TRUE, quote = FALSE, append = FALSE)
}

Thanks


text:
To find p-value of compounds (pairwise)
I have next SNPs data.frame structure
> str(SNPs)
'data.frame':   1703 obs. of  4 variables:
 $ group: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
 $ rs1  : Factor w/ 3 levels "D/D","I/D","I/I": 1 1 2 3 3 2 1 1 1 1 ...
 $ rs2  : Factor w/ 3 levels "a/a","a/b","b/b": 3 3 2 3 3 2 2 3 3 2 ...
 $ rs3  : Factor w/ 3 levels "G/G","G/T","T/T": 2 1 2 1 1 3 1 1 2 1 ...
 ...other rs

> head(SNPs)
  group rs1 rs2 rs3 ...other rs
1     A D/D b/b G/T
2     A D/D b/b G/G
3     A I/D a/b G/T
4     A I/I b/b G/G
5     A I/I b/b G/G
6     A I/D a/b T/T

For example, I noticed that rs5 4b/4a and rs6 G/G in group A very often occur together (see below). In group B they very often occur together too. So I want to know - is it statistical regularity or not.
I can create table with all pairs in both groups
> SNPs$rs5_rs6 <- paste(SNPs$rs5, SNPs$rs6)
> tmp <- table(SNPs$rs5_rs6, SNPs$group)
> tmp

              A   B
  4a/4a G/G   1  20
  4b/4a G/G  31  83
  4b/4a G/T  14  51
  4b/4a T/T   1   0
  4b/4b G/G  37 106
  4b/4b G/T  35 119
  4b/4b T/T  11  31

So, now I need compare (find p-value) group A and group B: 4a/4a G/G in group A vs 4a/4a G/G in group B, 4b/4a G/G in group A vs 4b/4a G/G in group B, 4b/4a G/T in group A vs 4b/4a G/T in group B and etc. Only this pairs/
How can I do that?
I need p-values for all
> tmp

              A   B
  4a/4a G/G   1  20 - p-value?
  4b/4a G/G  31  83 - p-value?
  4b/4a G/T  14  51 - p-value?
  4b/4a T/T   1   0 - p-value?
  4b/4b G/G  37 106 - p-value?
  4b/4b G/T  35 119 - p-value?
  4b/4b T/T  11  31 - p-value?

If you can restate your question in terms of explicit null and alternative hypotheses it will be much easier for people to help you determine an appropriate test to find p-values

@heathobrien chisq - it is my choice. Problem is not to choice test. Problem is to calculate it for pairs in tmp

If your question is how to apply a chisq test to the above table, the answer is chisq.test(tmp$A, tmp$B)

heathobrien, my question is how to apply chisq test to tmp and get 7 p-values:
compare4a/4a G/G in group A vs 4a/4a G/G in group B and get first p-value, compare 4b/4a G/G in group A vs 4b/4a G/G in group B and get the second p-value and etc.

OK. Now we're getting somewhere. See my answer

You can get standardised residuals from chisq.test and convert these to p-values, but it's not recommended. The frequencies of each genotype depend on the counts of all other genotypes, so the p-values aren't independent and they aren't corrected for multiple testing, but if you must:
pnorm(-abs(chisq.test(tmp[,c('A','B')])$stdres)))

hm... What other options are possible? whether the differences between the two groups were statistically significant for all pairs (rs5+rs6 and all pairs - 4a/4a G/G, 4b/4a G/G and etc)?

If I were you, I would just focus on the overall chisq p-value (corrected for the number of SNP pairs tested). The frequency of allele 2 depends on the frequency of allele 1 and the frequency of the genotypes depends on the allele frequencies, so I don't think it makes much sense to test each genotype independently

Note that all you're doing is looking at differences in linkage disequilibrium between groups. There are a lot of tools out there to that efficiently

chisq.test - not very good for values less then 5. I'v try to use fisher.test pnorm(-abs(fisher.test(tmp[,c('A','B')], simulate.p.value=TRUE)$stdres)), but I'v got the error - Error in abs(fisher.test(tmp[, c("A", "B")], simulate.p.value = TRUE)$stdres) :   non-numeric argument to mathematical function

o the p-values aren't independent and they aren't corrected for multiple testing, but if you must - I'm not must. I just think it is good solution to use chisq.test. I just need to know - is there any differencies in two grops

Fishers exact test doesn't give you residuals. Having values less than 5 is probably the least of your worries if you're trying to obtain p-values for each combination of genotypes.

In that case, I'm completely confused and do not understand what I can do with my data and how...

Let us continue this discussion in chat.


text:
What is the correct way of dealing with the analysis of data from flow cytometry?
I would like to detect a change in expression of a molecule present on a cell type by flow cytometry.
Assuming I am able to detect, using an antibody, a signal that represents the amount of the molecule I'm interested into. Also assume that the same experiment will have some noise (ideally near to zero) represented by an isotype antibody. What if the noise is somehow higher than a minimal threshold (i.e. = 1%)? How can I deal with the noise? How can I compare different samples if the noise is higher than the threshold?
I was thinking about removing those samples that have a noise above the 1% but is this ok, from a statistical point of view?
Another possibility is: can I subtract the noise from the signal and use this information for the final statistical analysis? Is this correct? Is it possible to subtract a signal from another? How is this done?
EDIT just to clarify a bit:
I have a pool of neutrophils which, in resting condition, do not express the molecule on their surface. When stimulated, these cells are supposed to show an increased amount of this molecule (again, on their surface). We would like to quantify the change in the amount of molecule (we used FACS).
We are able to quantify the quantity of surface-molecule but we also get some aspecific signal (the noise due to the isotype antiboby). What we try to do is to keep the isotype signal lower than 1% using a gate-based methodology. What if the isotype signal we get is higher than the threshold (i.e. >1%)? The antibody sometimes work some others don't. What should we do when it returns too much noise? Should we discard the data or should we perform a normalisation over the surface-molecule using the noise?

If your goal is to simply "to quantify the change in the amount of molecule," then you can do that at the population-average level by simply taking the average antibody signal for the different cell populations (stimulated vs. not stimulated). In that case, you don't have to throw out samples with noise above a threshold - since they are still part of the population and so should be used when calculating average.
If you want to quantify "the change in the amount of molecule" at the single-cell level, then things are a bit more complicated, and the best answer depends in part on things like the shapes of the distribution of signal from the stimulated and non-stimulated populations.
And related to that, when you say "noise above the 1%," do you mean that the Coefficient of Variation (CV) is greater than 1%? Or something else?


text:
Why is cutadapt discarding some paired end reads?
So, I am using cutadapt to remove two primers from some paired end reads.
I have 2 files for each sample, i.e. sample_R1.fastq and sample_R2.fastq.
I have primer_1 and primer_2 and both of them can be found in sample_R1.fastq and sample_R2.fastq. In other words, if I search for primer_1 and primer_2 in sample_R1.fastq, the sum of the matches is the same as the number of the reads in sample_R1.fastq. Same for R2.
To properly trim them what I do is the following, through R:
setOne <- paste("-g", primer_1, "-a", primer_1)
setTwo <- paste("-G", primer_2, "-A", primer_2)
system2(cutadapt, args = c(setOne, setTwo, "-n", 2,
                           "-o", outOne, "-p", outTwo,
                            inputOne, inputTwo, "--minimum-length", 1))

The command works fine except for one thing. On the one hand cutadapt should remove those sequences which, after the trimming, will have a length < 1, as specified by the command. However, some reads are removed even though they are longer than 1 and I would like to know why this happens, since I can't find a specific reference mentioning that cutadapt removes sequences.
If I am getting this straight, cutadapt only removes the primer sequence from the beginning of the read, then removes the equivalent characters in the quality field of the read that was trimmed and, finally, it returns the trimmed read with the quality field adjusted to the new length. If the length is lower than the threshold, the read is discarded. If it is discarded in R1, it will be discarded in R2 and viceversa.
I give an example of read that is removed, this one from R1:
@some:sample:number 1:N:0:1
ACCGTCACGTCATCAGCAGCATCGACTCCTGTTTGATCCCCGTGCCTGTACTCAAGTCAACCAGTACACGCTTTTCAAGCTGCCTTCTGCCGTATTTCACACTCAGCGCATTTGCAATTGGTGTTCTTTGCACCACATATCTATCAGTTAAGCCACCGCTACTCCGCCTACCTCGTACTGACTTAATCGACCGCCTACGCACCCTTTAAACCCAATTCGAGACCGTCAATTCAATCATACTTTAGTAAGTG
+
BBBBBBBBBBBBGGGGHEGHFHA@EGGGHHHHHHHHHGHHHHHHHHGGGGGHGHGHGGGHHHHHGDGGGFHGGGGHGGGGGGHGHGFHHHHHHHHHHGGGHHHGHGHHHHHHHGFFGGGGGGGGGGGGGGGGGFFBDFDFFFFFFFBFFFFFFFFFHGHGGHHHHHEGGGGGHHHHGHHHHHGHHHHHHHHHHHHGGGHGHHGGGGEFGGHHHHHHHHHHHHHGHHHHHHHHHHHFHFGHHHHGGGGGHHG

this one from R2:
@some:sample:number 2:N:0:1
CTGCTGCTATCACTAGCTAAAGCTTTAACGTAGGGCAATTGACGAAGGTGAAATGCATAGATATGTGATGACGGTGTAGCGTTGGCGCTGGGGTTTAAAGGGTGCGTAGTGACTTCACAGAACACCAATTGGCAGCTTACTAAAGCATGAGTACAGACGAGGTAGGCGGAATTGCGGTCGATTAAGTCAGTTGTGAAATACGGCAGCTTAACTGTTGGCAAGCGTTGTCCGGCACTGGTTTTATTCGAGGT
+
AAAAA@BFFFFFGFGCGGFDDC/CDHF2FGFHHG//<<CCC1FGHH1<1>11C<FF/0BFFFFFFF/:DEEEFFHGG0AEG/B9B./-.@9GHGEFC>EE<FHEHEGHDHGGFHGGBGGGHGGDGH0CDDGFFHHCC--.:CGCFG??GGAEGGHFHHBFFFEFB./B..;/;:/999=AAA../BB/;//////9/EFFCBDAA..990ADD=99/;/BFGGGGGGGGGFHHHHHHGHDGHHFFDHHHHG

and these are the primers that I am removing:

ACCGTCACGTCATCAGCAGCATCGAC
CTGCTGCTATCACTAGCTAAAGCTTTA

matching the reads I provided. Clearly, if I trim them manually, the reads length is way greater than 1 so, why are these reads removed from the final cutadapt output file?
EDIT:
I might have solved the issue, but I asked on the cutadapt github issues page, just to make sure.


text:
How likely is it to find primers sequences in already-trimmed reads?
So, I am analysing 18S amplicon data (fastq files) to be able to eventually investigate the taxa composition.
After removing primers (using cutadapt) from both R1 and R2 (assuming primers are found at the beginning of my sequences), I am still able to detect, apparently at random, a few reads that contain my primers (i.e. a part of the sequence matching the primer) in, again apparently random, positions.
I am able to detect almost all my primers (16 out of 24, both forward and reverse) in at least one sequence and I am wondering how common (or how weird) this is. I guess it is important to consider that my primers contains a few, i.e. from 1 to 3, wild-card bases that may match multiple other bases (apparently called whobble bases), which makes it easier to match a random sequence (there are multiple possibilities for an actual match).
Problem is that the presence of primers in the sequences may affect downstream analyses (e.g. with DADA2 they are very clear about the fact that primers have to be removed) in several ways and I was wondering if these sequences I am finding may be interfering somehow of if I can live with it.

The primers will be at the 3'-end of the reads. Make sure you check for them there.

@DevonRyan yes, that I know and that I did so no primers at 3'. I guess it just happens to find them inside the reads, as random matches, but I wasn't expecting that many actual matches. I guess it's statistics, then.

Interesting, have you tried to map a read with a sequence matching a primer to the reference, or even simply blasting it? Does the reference contain the "primer sequence"? ("" as we don't know if it's a homologous sequence or not)

@KamilSJaron apparently this is something that should not be happening and I'm trying to understand why it does.


text:
how to annotate the effect of a genomic variant to predict what the resulting protein is?
Spoiler: This is probably a theoretical question. But, of course, if you have any tool that does this I will be extremely happy.
So, I have WGS data for an organism O, both wild type (WT-O) and modified (M-O). I was able to assemble the reads using the O reference genome. Then, I was able to compare WT-O against M-O. I finally have the variants files (vcf), which I annotated with SnpEff and VEP.
Now what I need to do is to evaluate the effect of the variants at protein level. What I know from the vcf files is which nucleotides were affected, their position, which chromosome, gene, and whatnot (in most cases, although this is not always true for all the variants).
Question is: how to translate these variants from nucleotides to amino acids (AAs) and then, eventually, compare (BLAST?) the AAs sequences containing the variants to check whether they match with any known protein?
EDIT:
so the final goal would be to know whether the variants are causing some weird protein to be produced which may have a deleterious effect. It can be something like a toxin, or an allergen that are normally not present but that can appear "by random chance". I need to make sure that nothing like that happened and to be sure i need to know which proteins may result from the variants.

Doesn't VEP provide the effect at the protein level? At the very least it should give the reference and modified amino acids, right? Can you please show us a few lines of your VCF so we know what you have to work with? Also, I assume you don't have a budget for paid tools here, right? The company I work for has an API that can translate variants to HGVS protein notation and from that you can write a script that can build the protein sequence. But it's a paid service and you would still need to write your own script to get the full protein sequence.

Can you also give us some background on the biological question you are trying to answer? Why would you want to compare the proteins? It seem like a strange thing to do, I can't really think of why that would be useful. What is the final objective here?
What information do you want to get from this comparison? And what is this organism, does it have splicing or could you just translate the DNA directly?

@terdon, thank your for the link. i updated the question and i hope the final goal is more clear, now. also...i am not sure i can provide you with a sample from the VEP annotated file but i'll look into it.

@terdon, maybe i have what i need in this VEP file, but it is not very clear to me how it works. i am quite new to this field of variant calling. i am looking into it.

Just open the vcf file you have annotated with VEP, and show us a couple of variant lines with the annotations. Feel free to change the position/chromosome if you fear you cannot share it (although a couple of variants cannot be sensitive information). I want to see what kind of protein-level information you have. But the rest of your experiment doesn't really make sense to me. It sounds like you want to investigate the pathogenicity of the variants. There is no reason to try and re-create the protein and search for it, especially since the variants will often only change one amino acid.

@terdon, i understand but what if i need to have written proof of that? i need to make 100% sure that something unexpected happened.

You don't need them all, you just need a structure on PDB then the question is solvable.

@M__ can you explain what you mean?

You need a PDB file, without it I don't have a solution.