text:
What does liability mean in GWAS heritability?
I am reading about GWAS in heritability. They usually say to calculate heritability on a liability scale but after searching for this word "liability", I still don't understand clearly what does liability mean? Does anyone have any easy-to-understand explanation regarding this term? Thank you very much!

When you use linear mixed models to estimate heritability you assume that the underlying trait is normally distributed which is called a disease liability scale.
For continuous traits this is not a problem but for binary traits, this becomes an issue because you have a 0/1 value for a phenotype and usually there is a higher proportion of cases in the study sample than the general population prevalence of disease, which leads to an ascertainment bias.
The population prevalence also varies between populations, for example the prevalence of malaria in one continent is vastly different from another.
So to make the heritability estimates comparable and also to correct for the ascertainment bias, the observed scale heritabilities for dichotomous traits are transformed to the liability scale taking population prevalence and the case proportion into account.
Some good papers for reading:
Falconer D.S. The inheritance of liability to certain diseases, estimated from the incidence among relatives. Ann. Hum. Genet. 1965;29:51–71
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3059431/


text:
Database for proteome-wide predictions of protein structures
Accuracies of protein structure predictors have improved quite a lot in recent years. Algorithms such as Rosetta have gotten robust enough to predict structures of large number of proteins. However, I could't find any initiative to make a database for proteome-wide predictions of protein structures or may be integrate the predicted structures in existing databases (PDB or Uniprot).
If any such exists, please let me know.

So, your question is if there is a database of predicted structures? (What would be the utility of this database? How would it be up to date/updated?)

Yes.
There could be many uses of predicted structures of all the proteins in the proteome, given that the predictions are accurate enough. One could use that resource to make proteome wide inferences about protein structures, which we can not do now with structures of fewer than ~30% of the proteins available.
In terms of advantage, the users would not have to run the prediction model themselves, so it would save time and effort for them.

Regarding updation of the database, I would imagine that would a minor task if there is a prediction model available. Its just about rerunning the updated prediction model again.

I'm not sure that the accuracy is "good enough". Enough for what? Also the folding of a protein changes with the environment. If it seems so easy and no one have done it there might be some "hidden" problems on it. Perhaps the computational cost of those predictions is too expensive for any single (or a consortium) group to try.

There are three initiatives I know of to have a go at this:

PConsFam, which collects data from this paper.
The Baker group's metagenomic study which you mention in the question.
The recent DMPfold work from our lab. This compares its results to the above two studies and discusses the effect on various model organisms.

The second and third of these don't provide a database per se, but the models are available for download and are linked to the Pfam database.
EDIT: things have moved on significantly since I first wrote this answer. As other answers mention, the AlphaFold Protein Structure Database is now the best resource in this area.

DMPfold looks very promising.
Hoping that the work gets published soon so we could get to use the data.

Thanks, it's currently under review. Though all the data is already available, see https://github.com/psipred/DMPfold#data.

The paper is now published, see https://www.nature.com/articles/s41467-019-11994-0.

By far the most accurate model right now should be AlphaFold, and there are open structures across several proteomes available at https://alphafold.ebi.ac.uk/.

Updated
As of 2021, the EBI alphaFold2 is the best resource. There are some caveats, which I have outlined in this blog post.
But majorly, it contains only monomeric predictions and not oligomers and complexes. As a result, for protein with a close homologue with a solved structure, it may be more handy to use a threaded structure as found in the SwissModel database (see original answer below) because it may be oligomeric and/or one can copy over ligands or binding partners if deemed suitable.
One can make one's own predictions with AlphaFold2, who have released a handy notebook. Also one can make complexes and oligomers with some tweaks or with the heavily improved notebooks found in ColabFold GitHub repo.
Original Answer
It is worth saying that for model organisms, Expasy Swissmodel is the "classical" database of predicted proteome structures. It has several model organisms and gets updates regularly. The models are threaded only and the algorithm fell behind some time ago —their website however is modern and sleek. If using Rosetta afterwards, one should do a few cycles of FastRelax with a movemap preventing backbone movements or CA restrained, followed by a couple of cartesian bond/angle fixing cycles: the structure tends to blow up if energy minimised otherwise.
However, if one had a specific protein in mind predicting with RosettaCM, the I-Tasser server or Phyre server is by the recommended option. And actually, the Phyre server has a dataset of the precomputed human proteome.
About the PDB, the database has been purged of all in silico models and contains only empirical structures (NMR, X-ray, neutron diffraction and cryo-EM).

I am using DMPfold to predict protein structure of a 382 residue protein sequence. All goes wel, and even the files are produced but during the "seq2maps.csh", I get this error in between:
Running PSIPRED & SOLVPRED
Running PSI-BLAST with sequence c600m2-tail-before-lysin-sequence.temp.fasta ...
Segmentation fault (core dumped)
[blastpgp] WARNING: Unable to open BLOSUM62
[blastpgp] WARNING: BlastScoreBlkMatFill returned non-zero status
[blastpgp] WARNING: SetUpBlastSearch failed.
BLASTP 2.2.17 [Aug-26-2007]

Reference: Altschul, Stephen F., Thomas L. Madden, Alejandro A. Schaffer,
Jinghui Zhang, Zheng Zhang, Webb Miller, and David J. Lipman (1997),
"Gapped BLAST and PSI-BLAST: a new generation of protein database search
programs",  Nucleic Acids Res. 25:3389-3402.

Reference for composition-based statistics:
Schaffer, Alejandro A., L. Aravind, Thomas L. Madden,
Sergei Shavirin, John L. Spouge, Yuri I. Wolf,
Eugene V. Koonin, and Stephen F. Altschul (2001),
"Improving the accuracy of PSI-BLAST protein database searches with
composition-based statistics and other refinements",  Nucleic Acids Res. 29:2994-3005.

Query= C600M2_00100 tail protein [Salmonella phage BPS17S6]
         (377 letters)

Database: c600m2-tail-before-lysin-sequence.a3m
           46,978 sequences; 11,323,713 total letters

SearchingFATAL: Error whilst running blastpgp - script terminated!
Running PSICOV
Running FreeContact
Running CCMpred
 _____ _____ _____               _
|     |     |     |___ ___ ___ _| |
|   --|   --| | | | . |  _| -_| . |
|_____|_____|_|_|_|  _|_| |___|___|
                  |_|

using CPU (10 thread(s))

Reweighted 30404 sequences with threshold 0.8 to Beff=22917.9 weight mean=0.753779, min=0.0232558, max=1

Will optimize 62686429 32-bit variables

iter    eval    f(x)        ║x║         ║g║         step
1       1       1.50338e+07 152357      3.021629e+11    1.55e-06
2       1       1.47369e+07 152356      2.1447226e+11   1.18e-06
3       1       1.44292e+07 152355      1.4646483e+11   1.22e-06
4       1       1.41268e+07 152353      1.0388224e+11   1.5e-06
5       1       1.38295e+07 152358      7.337667e+10    1.92e-06

Can you please let me know if I need to rerun seq2maps.csh because of this error? And, how this error can be solved? Thank you!

I will answer where you have posted this as a separate question: https://bioinformatics.stackexchange.com/questions/9345/dmpfold-error-while-running-protein-structure-prediction

source: Stack Exchange

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

sort: votes


text:
How to install genometools with anaconda?
I tied to install genometools with anaconda:
conda install -c bioconda genometools

Though, after running that successfully the command gt is not available in my conda environment.

The program 'gt' is currently not installed. You can install it by
  typing: sudo apt install genometools

While conda tells me the program is installed. Are there multiple packages called genometools? What is the difference and can I install the package providing gt with anaconda at all? There is also a package called genometools-genometools, but it only has two downloads. So, I am not sure that is what I need.

Possible duplicate of Difference between genometools installed with anaconda vs apt-get?

genometools-genometools seem to be the current version. By installing the genometools package I received the same error like you. I guess it's because the binary isn't copy to bin during installation.
genometools-genometools seems to work fine. And it has 2k (2000) and not 2 :)

Maybe I was looking at the wrong channel. Thanks!

genometools-genometools worked for me, while 2 other options (genometools and -c leylabmpi genometools) did not.

Hi thanks Sara could you clarify your installation code please?  it appears you are designating channnels within Anaconda

Hi Michael, the example with the channel specified was a "fix" that was suggested elsewhere. What worked for me was "conda install -n GT genometools-genometools" where I specified the environment I wanted to install in (GT) but not the channel to use.

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


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:
Output the linear predictor from a stratified cox model?
I have fitted a cox model on a pooled dataset of multiple studies, say studies A, B, and C. As these studies have different baseline hazard functions, I stratified for 'study' in the model.
Like this:
df$SurvObj <- with(df, Surv(event_rd, event == 1))

fit <- coxph(SurvObj ~ cov1 + cov2 + cov3 + cov4 + strata(study),data=df)

Now, I want to assess whether I can use the beta coefficients from the above model to predict the event probability in a new study, study D.
The baseline hazard in study D is different from those in studies A, B, and C.
When I have
predict(fit,type="lp",newdata=studyD,reference="strata")

I get this error:

Error in model.frame.default(data = studyD, formula = ~hba1c + sbp +
: factor strata(study) has new levels study=4

Why does R require that I match the strata of study D with those in study A, B, and/or C when I only want to output the linear predictor (type="lp")? As the linear predictor is independent of the baseline hazard function in this case.
I suppose it should be possible to extract the linear predictor for individuals in study D and then manually calculate the event probability using the baseline hazard of study D.
Does anybody know how to do this?
I could not find any question like this so maybe I am thinking the wrong way.

I would recommend trying out Stats.SE (aka CrossValidated) as some folks there might be of more help.

Unlike wikipedia, the linear predictor (LP) in coxph is made dependent from the baseline hazard function $\lambda$. For example, if a model has two predictors x1, x2 stratified by sex, the linear predictor would be like:
$$LP = log (  \lambda _{sex}) + \beta _{1}x _{1}+\beta _{2}x _{2}$$
Where the $log(\lambda _{sex})$ can be obtained by predict() with x1=x2=0.
We could calculate the LP without strata information (sex) by hand using the $\beta _{1} \, \beta _{2}$ provided by coef(fit) and a customised $\lambda _{unkown \, sex}$.
Below is a simulation comparing the results from predict() and by hand coef(fit).
## Toy data
df <- list(time=c(4,3,1,1,2,2,3),
          status=c(1,1,1,0,1,1,0),
          x1=c(0,2,1,1,1,0,0),
          x2=1:7,
          sex=c(0,0,0,0,1,1,1))

## Fit a stratified model
fit <- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df)
coef(fit)
#         x1          x2
# 0.79451881 -0.01917633

## Baseline linear predictor
lambda_sex1 <- predict(fit, newdata=list(x1=0, x2=0, sex=0))
# -0.746578
lambda_sex2 <- predict(fit, newdata=list(x1=0, x2=0, sex=1))
# -0.1497816

## by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=0))
# 0.009588167
predict(fit, newdata=list(x1=1, x2=2, sex=1))
# 0.6063845

## by hand
lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.009588167
lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.6063845

Finally, to calculate "a" linear predictor for an unknown sex:
# prevented by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=2))
#Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  :
# factor strata(sex) has new level sex=2

# by hand, assuming lambda for sex 2 is 0
0 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.7561661


text:
Plotting coverage of annotation over collection of region
I'm trying to plot "meta" coverage of annotation: i.e. features (eg. gene class) over certain regions. It is similar to read coverage plots over gene body, except my input is two bed files (both in BED6 format) - (A) one containing the regions for which the meta plot is needed, and another being the list of genes or any feature whose coverage over the regions in (A) is needed.
Is there any package or tool which can create such plots (my domain is limited to python but I can try to work with R) ?
Something akin to this but for whole gene body[1]:

[1]- Kelley, D, Rinn, J (2012). Transposable elements reveal a stem cell-specific class of long noncoding RNAs. Genome Biol., 13, 11:R107.
Edit: The BED file containing list of regions looks like this:
chr3    39218734    40053659    region1 0.92426187419769    +
chr4    140163762   140453127   region2 0.896103896103896   -
chr7    40549151    41205036    region3 0.986072423398329   +
chr8    81291743    81963246    region4 0.94184168012924    -
chr9    12284032    12539789    region5 0.95539033457249    -

And the bed file containing features to be plotted looks like this:
chr3    39218100    40053200    LINE    1   +
chr4    140163962   140453027   LINE 1  -
chr7    40549002    41204999    SINE 1  +
chr8    81291143    81963846    LTR 1   -
chr9    12284332    12539720    LTR 1   -

I am not familiar with this field, so I would like to know the terms I should avoid/are not good enough. What kind of keywords have you searched?

@llrs, I'm not exactly sure what you mean by "terms I should avoid/are not good enough." - but I believe coverage plots are one class of visualisation which are more similar to binned histograms in general. Except a google search would yield in methods related to plotting coverage of transcripts or reads acquired from sequencing experiments [keywords: coverage plots, bedfiles, metaplots]. The key difference here is that I'm interested in scaled coverage over whole gene body (TSS-TES) instead of positions relative to TSS - similar to deeptools plotProfile.

If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:

Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
Convert those to bedGraph (e.g., awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed > LINE.bedGraph).
Use UCSC tools to convert those bedGraph files to BigWig.
Install deepTools and run computeMatrix reference-point -b 2500 -a 20 -S LINE.bigWig SINE.bigWig LTR.bigWig -R Regions_Of_Interest.BED -o foo.mat.gz
Make a profile plot with plotProfile (plotProfile -m foo.mat.gz -o foo.png --perGroup)

If you do have overlapping regions you'll need to first make a disjoint set of intervals (you can use bedops for this), find the coverage of them (e.g., bedtools coverage) and then continue on with that.

I still have to try it, but this might just work. I will mark it as answer unless something goes wrong. I would use scale-regions instead of reference-point since I'm interested in the whole region length.

There is a reasonably nice looking tutorial on metaplots in R here: https://rpubs.com/achitsaz/94710

Thanks for answering this, usually we'd like to have the answer self-contained. Perhaps you could improve the question by including the relevant packages/code from that website into your answer.


text:
BPGA - functional analysis
I have been using BPGA to perform pan-genome analysis on my set of isolates. Now, I need the functional analysis as well. However I am having problem with the function, as the program says I am missing two files: usearch.exe (which I am sure is already present) and cog.udb, which I do not know what it is and where to find it.
Any suggestions on how to solve the error?
Thanks

Welcome to the site Rob. Could you provide more information about your settings: You seem to use Windows, could you copy and paste the error message you get? Did you follow the user guide installation ?

The error is "File/s: cog.udb & usearch.exe not found. Copy missing files to current folder."
And yes, I followed the user guide for installation. Thanks!

Did you follow the installation steps of the linked user guide?

Yes, I followed the instruction


text:
Which software is used to run Molecular Dynamics simulation of a small RNA hairpin?
I am currently working in a project that needs MD simulation of a small RNA hairpin. I am currently trying it on GROMACS but it seems to be protein based, hence am facing several issues.
Is there any software that is designed for the same?

Welcome to the site @AKB (I hope you register here also). How long is your small RNA ? How long do you need your simulation to go on? Also if you could provide how have you searched for those programs we could help you to refine the search query too

Google found this review that might be of interest: https://pubs.acs.org/doi/10.1021/acs.chemrev.7b00427
Here is a link to a free full-text version:
https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/29297679/
The authors "here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA–ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field."
Hope this helps.

I've added a link to a free version, and provided an abstract from the abstract.

AMBER has force fields specifically for viewing RNA free energy and structural dynamics. Specifically, AMBERff99 is useful for RNA analysis. It is fairly easy to use and a comprehensive software for your analysis. It has plenty of tutorials from basic to advanced usage, as well as many publications. A specific publication related to your question can be found at
doi: 10.1016/j.bpj.2009.10.040
But in general, the J Phys Chem B is relevant to MD studies.


text:
Output from vcftools missingness
I'm new to data filtering on vcf data and vcftools.
I performed variant calling on my dataset, CHR22, homo sapiens. I'd like to remove sites that are missing in more than 5% of individuals.
vcftools --missing-site --vcf updated_ids68.vcf

This gives me a file with out.lmiss. That is missingness per locus.
I had 68 individuals in my input vcf file.
CHR     POS     N_DATA  N_GENOTYPE_FILTERED     N_MISS  F_MISS
22      16848278        106     0       30      0.283019
22      16848492        68      0       68      1
22      16849180        69      0       67      0.971014
22      16849229        68      0       68      1
22      16849376        133     0       3       0.0225564
22      16849476        132     0       6       0.0454545
22      16851734        126     0       36      0.285714
22      16852588        123     0       13      0.105691

I'm unable to understand the output here. What does third column - N_DATA tell?
133, 126 doesn't tell that this is number of individuals if my input had 68 individuals only.
VCFtools (0.1.15)
fileformat=VCFv4.1

Humans are diploid, so you can expect to see up to 2*N (2*68=136) alleles, so N_DATA is the number of observances of that allele.

Furthermore, it also looks like N_MISS + N_DATA = 2*68

My initial answer mentioned that, but that doesn't seem to hold true for all rows. Consider the following rows from their table:
22      16849476        132     0       6       0.0454545
22      16851734        126     0       36      0.285714
The N_DATA and N_MISS instead sum to 138 and 162, respectively.

Right. That is surprising. What could cause that?

@spvensko Great catch.  My data are structural variants, maybe that is why?  Thanks.


text:
Pathway plot in GeneTrail data base
So I'm doing a pathway analysis using GeneTrail database. One thing I am not able to get is this similarity measure being calculated and plotted:

I want to know how this is being made, what is the output being plotted because each gene is may be involved in multiple pathways.
I want to be able to make it in R again. Any suggestion or help would be appreciated. How to plot the enrichment result in this form ?

This is a a symmetric matrix with the pathways in both axes. So for each pair a similarity index is calculated (usually it takes into account the number of genes in each pathway and the total number of genes in both pathways). There are several similarity indexes available. A typical one is Jaccard or Sørensen-Dice index, usually 0 means no similarity and 1 equal pathways.
You can do this kind of analysis in R using BioCor (Disclosure: I am the author of this package) which uses the Dice similarity. You will need the name of the genes of each pathway, after calculating the similarity of pathways using mpathSim, you can then plot it with a heatmap (using the heatmap function).

yes i remember  from my last post you have given the link im going to use it and let you know ..i think this looks bit more comprehensive to put all the pathway together..if i have any issue with your library i will ask you..


text:
Iupred definition of long/short form disorder
I need to predict disorder of proteins and here there is a description of Iupred in its two variants, long and short-disorder. This is a manuscript defining the method. I couldn't find a precise definition of short or long-disorder.
I need to predict disorder in 30-mers, but I expect that as with secondary structure prediction, the longer the input sequence the more accurate the result.
I will appreciate if anyone has experience with this software and can share ideas or suggestion on how to proceed.


text:
Restricting match output for multiple sequence alignment using mafft?
So I aligned roughly 5k sequences and I got my output using mafft. However, I want to restrict the output to only present conserved segments of a specific length (let's say 25bp). Does mafft have parameters to set this that I seem to be missing?

Q - Is there some parameter for this?
No, there is no parameter for this. Getting the conserved regions is something that you will have to perform post-alignment.
How to get the conserved regions?
You will need to define what exactly you mean by conservation and calculate a conservation score for each position of the alignment.
For instance you can look at each column of the alignment and calculate the percent of queries that had the same base/AA at that position. If using AAs you can define conservation more loosely depending on what properties are important to you, or using a scoring matrix.
Once you have a co-occurrence matrix you can define the regions of interest based on a conservation threshold that you define.
Is there software to find conserved columns?
There is a program MstatX which outputs conservation scores of individual columns.
How MstatX works
The example of MstatX provides this example alignment. Each column is progressively less conserved than the last. The last column is conserved, but has some gaps -- this lowers the conservation score.
   col          1
 seq1  1234567890
 seq2  DDDDDDIPDL
 seq3  DDDDDDIPVL
 seq4  DDDDDDIPYL
 seq5  DDDDDDIPA-
 seq6  DDDDDDLWT-
 seq7  DDEDEELWK-
 seq8  DDEDEELWP-
 seq9  DDEDEFVSR-
 seq10 DEEFFFVSH-

The conservation scores for each column were:
1       2.71272e-08
2       0.243433
3       0.312789
4       0.243433
5       0.483364
6       0.480391
7       0.498757
8       0.498757
9       0.973514
10      0.267319

How to get conserved regions
The python program below will output all windows of conserved regions when using a fasta alignment and the default output from MstatX.

Positional args

1. MstatX_alignment_example.mali, is the fasta alignment file.
2. output.txt is the default output of mstatx when called with mstatx -i MstatX_alignment_example.mali.
3. 0.4 is the minimum conservation score
4. 3 is the minimum number of bases that must have that conservation score.

Calling this program with the following command on the above example:
python print_conserved_regions.py MstatX_alignment_example.mali output.txt 0.4 3 prints:
{1,4}

Meaning that there was only one window at least 3 bases long in which all positions had a conservation score below 0.4. This window was positions 1 through 4 in the alignment.
#call this file cons_windows.py
import sys
from Bio import AlignIO

align_filename = sys.argv[1]
cons_filename = sys.argv[2]
cutoff = float(sys.argv[3])
minchar = int(sys.argv[4])

def parse_conservation(filename):
    cons = {}
    with open(filename, "r") as f:
        for line in f:
            if line.strip():
                temp = line.split()
                cons[int(temp[0])] = float(temp[1])
    return cons

align = AlignIO.read(align_filename, "fasta")
cons = parse_conservation(sys.argv[2])

i = 1
while i < max(cons):
    if cons[i] <= cutoff:
        j = i
        done = False
        while j < max(cons) and not done:
            if (cons[j] > cutoff) and (j-i >= minchar):
                print("{{{},{}}}".format(i,j-1))
                done = True
                i = j
            j+=1
    i+=1


text:
How I know which gene is a good predictor in this neural network or not?
I have 25 highly differentially expressed genes among and patients to chemotherapy.
I have made a neural network of these genes. The accuracy of the model is 0.73 but I don't know from 25 genes how I could select the best predictors or remove the weaker ones.
> # Convert column from /to 1/0
> = as.numeric(data$)-1
> data1 = cbind(,df)

> head(data1[,1:4])
   ALB CALML5 CALML6
A2       0 4.64   5.77   3.16
A3       0 7.12   9.13   5.13
A4       0 7.84   7.09   4.85
A6       0 7.90   9.64   6.52
A7       0 5.84   6.26   0.00
A8       0 4.83   5.77   2.15
>
> library(caTools)
> set.seed(101)
> # Create Split (any column is fine)
> split = sample.split(data1$, SplitRatio = 0.70)
> # Split based off of split Boolean Vector
> train = subset(data1, split == TRUE)
> test = subset(data1, split == FALSE)
> feats <- names(data[,-1])
> # Concatenate strings
> f <- paste(feats,collapse=' + ')
> f <- paste('TG ~',f)
> # Convert to formula
> f <- as.formula(f)
> f
~ ALB + CALML5 + CALML6 + CDK6 + CHGA + CREB3L3 + CSF3 +
    CXCL5 + CXCL6 + CXCL8 + DKK1 + EGLN2 + FGF20 + FGF4 + IL11 +
    IL1B + IL6 + ISG15 + KLK5 + KRT17 + MMP3 + MMP7 + PKP1 +
    PTGS2 + S100A7A
> #install.packages('neuralnet')
> library(neuralnet)
Warning message:
package ‘neuralnet’ was built under R version 3.5.2
> nn <- neuralnet(f,train,hidden=c(25,2,1),linear.output=FALSE)
> View(test)
> # Compute Predictions off Test Set
> predicted.nn.values <- compute(nn,test[2:26])
> # Check out net.result
> print(head(predicted.nn.values$net.result))
             [,1]
A4  0.26890819169
B5  0.94836576493
B7  0.19753554824
B8  0.99790838807
B9  0.99764782073
B12 0.05230351043
> predicted.nn.values$net.result <- sapply(predicted.nn.values$net.result,round,digits=0)
> table(test$Private,predicted.nn.values$net.result)

    0 1
  0 4 3
  1 1 7

> accuracy<-(table(test$Private,predicted.nn.values$net.result)[1,1]+table(test$Private,predicted.nn.values$net.result)[2,2])/(sum(table(test$Private,predicted.nn.values$net.result)))

> accuracy
[1] 0.7333333333
>

Usually what I've heard is that one does it so iteratively, so you have now 25, you remove one and see what happens and then remove another one and so on...

Thank you, but what happen to what? I mean in regression I am checking pvalue of coefficient of variation but here I am not which predictor is more influential

Well, you are evaluating the accuracy, so this is your metric you want to follow/to maximize (although I wouldn't forget about specificity).

Yes, there is a formal solution. You want a heat map of the layers weighting and this is the information sought. It is more difficult to interpret than other supervised learning algorithms, but that is the nature of neural networks.
I know Sci-test Learn (Python), but I don't know R for neural networks, so I can't advise you on the code, only that it will be a few lines of code and a graphical (relative) output. I imagine Tensorflow will have a really brilliant heat map layer weighting.


text:
How to correctly call a VCF file using damaged DNA? (IonTorrent & FFPE)
EDIT: I am updating this question to make it more specific to my issue.
For context - original question prior to edit:
How do I obtain a deamination metric when doing the variant calling using the IonTorrent variant caller, and secondly, how do I correct my called variants for deamination to ensure that these don't provide false positives for the downstream analysis?
Background
I'm calling variants using the IonTorrent TorrentSuite on DNA which has been sequenced from formalin-fixed paraffin-embedded (FFPE) tissue. This has a major issue in that without addition of uracil-N-glycosylase, some of the Ts in the original DNA are deaminated to Cs, which upon sequencing and calling variants can show up as mutations, either as T>C transitions or G>A (from the opposite strand, due to PCR in the library prep). I do not have any idea how long these samples were stored without UNG before sequencing.
TVC (Torrent Variant Caller) gives a deamination metric (essentially, sum of T>C and G>A variants over all variants called), and for our samples, the highest value seen is ~0.92. Naively postprocessing the variants show that for these samples, C>T/T>C transitions overwhelm the remaining variants among my samples.
Application
My use case is this: these are medical samples, which have been inspected by a pathologist (hence the FFPE treatment), and I want to determine which variants are predictive of outcome, hence I have two potentially contradictory goals: reduce false positives and capture the rarer variants which may hold predictive power.

Question

Given the IonTorrent variant calling pipeline (sequencing > BAM file > TVC > VCF file with deamination statistic):

(Q1) Is there a way of correcting the output VCF to remove transition and deamination errors?
(Q2) Alternatively, is there a set of filters to use in bcftools to reduce the effect of the errors on how variants are called?

what is your intended downstream use of the reads? It may not be necessary to correct them at all.

I want to take my samples and make inferences about which variants correlate with differing responses to treatment in a more clinical setting, so not having any signal swamped by this C>T T>C noise is what I require

Sites that were deaminated are often called abasic sites.
Q1: How to obtain DNA damage metrics?
You might want to check out mapDamage - this software is used to perform these sorts of analyses and produce plots of types of DNA damage along sequencing reads.
Q2: How to correct for abasic sites and other DNA damage errors?
You can use your choice of kmer-based or alignment based error correction software. The errors will be random and if the sequencing depth is sufficient the erroneous kmers at the ends of reads will be corrected. However, depending on what you are trying to do downstream with your reads you may be better served by not correcting the reads, then dealing with transitions/damage later.

That's fine but I'm trying to combine some older IonTorrent called variants with variants I've had to call myself, so I don't have access to the original mapped read files for all the samples. I was essentially asking for which settings to try to get the deamination estimate, and how to correct for it once the vcf's have been produced. I want to take my c. 151 samples and make inferences about which variants correlate with differing responses to treatment in a more clinical setting, so not having any signal swamped by this C>T T>C noise is what I require

I have since found the deamination statistics as provided by TVC (the TorrentSuite variant caller) but still would like to know how to deal with the number that is given and how to filter my VCFs to ignore the false positives occurring due to this mechanism.

Having inspected one of my samples (with estimated deamination metric, the proportion of called SNPs consisting of C>T or G>A transitions, of 0.837), this looks like a straightforward filtering procedure. Looking at the minor allele fraction (AF) and the quality score (-10 log10(p-value)) per type of SNP:
it can be seen that the C>T and G>A transitions likely caused by deamination (and PCRing up the deaminated DNA to get the libraries) all cluster in the bottom right of the plot - low AF and low QUAL, so setting thresholds to the left (higher QUAL) and above (higher AF) the groups will remove these from the called variants.


text:
How to score how densely bases are distributed along a given oligo?
Suppose I have a short chain string of oligos, and I want to rate them from 0 to 1 based on how clustered the "A"s are in respect to the chain. For example:
"AAAAACCCCC" #should score as 1 for being as clustered as it possibly could
"ACACACACAC" #should be scored as 0 for being as thinly distributed as it could

Current idea is to strip both ends of A's, then remove first and last characters, then
gaps = ''.join(str(x) for x in deck).split("A")
listofCs = [x for x in gaps if x]

to join it all together as a list of how big the gaps are (i.e. the "C"s), and simply count the average size
average = sum(map(len,listofCs))/len(listofCs)

But I feel like there should be a much more elegant solution to this. It feels messy. How can I improve the scoring?

I would look at how much consecutive the different elements are in the sequence. But how would you rate a sequence like AACCAACACCAACC ? Is it for DNA or protein sequences (how much elements do you expect, 4, [5 with an N], 21) ?

You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As.  You might also need to scale by the sequence length, if you have different length sequences that you want to compare.
Edit: also, there might be a correction to make if one or both ends of the sequence are an A

Interestingly idea. Can you think of any edge cases that would cause problems for this scheme?

Zero As would throw a spanner in the works.

To define the base density score, or clustering score, of a sequence $S$ we can think of each position in the sequence as being part of a Markov process. For the sake of symmetry around a single base we consider odd-order Markov models from 3-onward where we ask what is the probability that we see the character at $S_i$ in the $j$ bases from $[S_{i-j}, S_{i+j}]$. For each column we have length(S) -1, $|S|-1$, possible probabilities to calculate per window.
We can't look left or right of the edges of the sequences so the probabilities are skewed, but we can get a rough idea of what is going on.

For each base position we calculate a clustering score of $ \frac{1}{w}\Sigma_{j=1}^{j=w} -log_2 \frac{m(S_{i-j}, S_{i+j})-1 + \frac{1}{|S|}}{|S_{\{i-j, i+j\}}| -1 }$ where:

$w$ is the window size we can look at. The max is $|S|-1$ and the min is $1$. This just means that we can calculate a score in the context of the whole sequence or we maybe only care within some smaller window, like 5 bp flanking either side of the base in question.
The function $m(S_{i-j}, S_{i+j})$ is just the number of counts of matching bases in the window $j$ on either side of the base $S_i$. For example in the sequence gatAgga, if $S_i$ is the uppercase A and $j$ is 3, the function $m(S_{i-3}, S_{i+3})$ returns 3, since there is the A base at $S_i$ and the two additional a bases at $i-2$ and at $i+3$. We have to subtract one from this number because we don't want to factor in the matching base at $S_i$.
We then add $\frac{1}{|S|}$ to that number. This is a pseudocount that we need to add to prevent from taking the log of zero. This number is pretty arbitrary but I chose $\frac{1}{|S|}$ so that non-clusting regions' scores aren't too high.
We then divide that numerator by the length of $S$ from the little window that we just looked at minus one, $|S_{\{i-j, i+j\}}| -1 $. For example in the sequence gatAgga from above we only want to consider the bases gat gga. This is 6 bases.
We then take the negative log to convert the proportional scale to a log scale.
The summation just says that we do this for every symmetrical window size around the base in question then sum them up. For example we would calculate the ratio for tAg, atAgg, and gatAgga.
The $\frac{1}{w}$ just gives the average log clustering score of all the windows that we looked at. One nice thing about this is that the impact of each increasing window size has less effect on the ratio, so there is a sort of convergence that happens.

Implemented in python:
from math import log as mlog
s = "GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA"
assert len(s) > 0
pseudocount = 1/len(s)
#maxdepth = 5
maxdepth = len(s)
score = []
for i in range(len(s)): # iterate positions
    temp = []
    for j in range(1, maxdepth): # iterate depth
        rmin = max(i-j, 0)
        rmax = i+j+1
        num_matches_min_one = sum([1 for k in s[rmin:rmax] if k == s[i]])-1
        len_s_min_one = len(s[rmin:rmax]) -1
        sim = abs(-1 * mlog((num_matches_min_one/len_s_min_one) + pseudocount, 2))
        #print("i={}, j={}, rmin={}, rmax={}, sim={}".format(i, j, rmin, rmax, sim))
        temp.append(sim)
    print("{0:d}\t{1}\t".format(i, s[i]), end='\t')
    score.append(sum(temp)/len(temp))
    print("{0:.4f}".format(score[-1]))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
x = [i+1 for i in range(len(s))]
plt.bar(x, score)
plt.xticks(x, [k for k in s])
ax.set_yscale("log", nonposy='clip')
plt.xlabel('sequence')
plt.ylabel('-log2(cluster-ness)')
plt.title("cluster-ness score")
plt.savefig('output.png', dpi=150)

This calculates the matrix, prints out a tab-delimited file of the position, base, and log clusting score, then saves a plot of the results.
For the test sequence GGATAGGACATTTTTTTTTTGAGCTTTTATTTTCCGA and using the whole sequence as context, I got the following plot:

As you can see, the log clustering score is low around the things that are cluster-y and high around things that aren't cluster-y.
Note: OP requested that the score be scaled from 0 to 1. I instead used the log scale to make visualization easier.
Warning 1: If you plan on doing this for large sequences, it is really important to limit the maxdepth, or $w$, parameter. The time complexity for this algorithm is $O(|S|*w)$, so the number of calculations can blow up to $O(|S|^2)$.
Warning 2: It is also important to limit the $w$ parameter because the edge effects propagate through the entire sequence when $w=|S|$. Look at these examples built off of OP's initial suggestions of AAAAACCCCC and ACACACACAC. The clustering score never stabilizes when $w=|S|$, however it does stabilize away from the edges when $w=5$.
For 50 characters of ACACAC...:

For 50 characters of [25xA][25xC]:

The Markov chain approach is very interesting. You could deploy a maximum likelihood approach where instead of scoring the genetic distance within a homologue, the algorithm looks at its physically linked neighbour. The genetic distance of the transition point i.e. AAAA then CCCC is deducted from the final genetic distance so that AAAAACCCCCC would be zero. It would be fairly easy to script and use PAML or equivalent where the default phylogeny would be a star tree.
The challenge is to establish a Deep Learning approach. Its not impossible, you just assign pre-determined classes of oligos which fall into the definition you want and use these as your training set.


text:
Gene Ranking - signal to noise ratio used in GSEA-P algorithm?
I'm looking at Broad Institute's orignal GSEA-P algorithm R script which I downloaded here: http://software.broadinstitute.org/gsea/downloads.jsp.
I'm trying to adapt their GSEA.1.0.R script to process datasets that have 1 gene expression profile as for class 1, and N gene expression profiles for class 2. The bottleneck seems to be the GSEA.MetaboliteRanking() method. I looked and it appears they are using "signal to noise ratio" (SNR) metric to rank genes in patient profiles.
I was wondering how the signal to noise is calculated from gene expression profiles?
Note: I am using/modifying their script because it's the only implementation of the GSEA algorithm that I could find where they expose their code.
i.e., what is the mathematical reasoning behind the following lines of code? Note, there are little to no comments in the code, but the following variables, from what I can tell, are:
A: Matrix of gene expression values (rows are genes, columns are samples)
class.labels: Phenotype of class distinction of interest. A vector of
              binary labels having first the 1s and then the 0s

P = class.labels1 * subset.mask
n1 = sum(P[,1])
M1 = A %*% P
M1 = M1/n1
A2 = A*A
S1 = A2 %*% P
S1 = S1/n1 - M1*M1
S1 = sqrt(abs((n1/(n1-1)) * S1))

P = class.labels2 * subset.mask
n2 = sum(P[,1])
M2 = A %*% P
M2 = M2/n2
A2 = A*A
S2 = A2 %*% P
S2 = S2/n2 - M2*M2
S2 = sqrt(abs((n2/(n2-1)) * S2))

The noise is comming from the subset.mask, which is created above in a loop with the number of permutations.
for (r in 1:nperm) { #L90
....
        subset.mask[, r] <- as.numeric(c(subset.class1, subset.class2)) # L107

So by multiplying the random selection of subsets by the expression we get the "noise". Later on line 251 we get the ratio of noise/signal by dividing the signal from the contrast and the signal by the noise.
Those transformation on the middle seem to come from normalizing the signal in order to be able to compare the resampling, but I don't fully understand why they do this. Also note that the R scripts are outdated and no longer maintained, so they can contain errors.


text:
What reasons are there to choose Illumina if PacBio provides longer and better reads?
PacBio provides longer read length than Illumina's short-length reads. Longer reads offer better opportunity for genome assembly, structural variant calling. It is not worse than short reads for calling SNP/indels, quantifying transcripts. Sounds like PacBio can do whatever Illumina platform can offer.
Would that be any reason for going with Illumina if PacBio's long reads can do everything and more?

There are so many reasons why one might want to prefer Illumina over PacBio (also note that it's a false dichotomy, at least Oxford Nanopore is a competitive sequencing platform):

The first (IMHO and the most common reason) is still the cost of both sequencing and the instruments. Illumina can sequence a Gbp of data for \$7 - \$93. PacBio sequencing is according the same webpage \$115 per Gbp, however at our sequencing center it's ~$200. Though ONT might have already a cheaper solution. Edit, I just found a google sheat with prices that seems to be frequently updated, seems that the ratio still holds Illumina short reads ~10x cheaper than PacBio.
RNA-seq (i.e. analysis of a gene expression) is not possible with PacBio due to preferential sequencing of smaller fragments; shorter genes would always be shown to be more expressed. To be clear, it's possible to sequence RNA with PacBio (the keyword is iso-seq), but the analysis of gene expression is problematic.
It's way easier to extract fragmented DNA (concerns small non-model organisms; although recently a single mosquito was sequenced, so we can expect a further improvement)
other sequencing techniques as RAD-seq that allow genotyping with very little effort and cost, I have never seen anybody even considering using long reads for such genotyping
Genome profiling (assembly-less genome evaluation) based on kmer spectra analysis is not possible with PacBio data due to higher error rates. Conflict of interest: I am a developer of one of the tools for genome profiling (smudgeplot)
DNA in Formalin-Fixed Paraffin-Embedded samples is fragmented anyway (~300bp), therefore, there is no point in sequencing them with more expensive long-read technology (contributed by @mRoten)
I bet there will be a plenty other applications I am not aware of

Please add that FFPE sample DNA is fragmented (~300 nt), so Illumina short read sequencing is still has the lowest cost per base than PacBio. Even ConcatSeq (https://www.nature.com/articles/s41598-017-05503-w) only increase the read count by a factor of 5. Illumina is still the winner.

I turned the answer to community wiki so it's easier for people to add/updates reasons for Illumina. As this is not complete not a stable-in-time answer.

Many analyses performed on Illumina machines these days require large numbers of reads. For example, most analyses in ChIP-seq, RNA-seq, ATAC-seq etc, need 10s or even 100s of millions of reads for the statistics to work out properly.
But this isn't limited to just sequencing as an assay experiments. High depths are important for things like somatic variant calling as well. Or simply sequencing of medical gene pannels, where you might only have 20 amplicons of a 500bp each, but need to do it on 100s of patients as cheaply as possible.

Three reasons for Illumina:

Much better for a large number of samples (easily handle 96 samples).
SNP calling is much better - much greater depth
Hardware costs, an Illumina MiSeq machine is cheap

PacBio SNP calling has been improved towards the standard of Illumina by DNA modification to the inputs to be sequenced but it depends on the application.

As others have mentioned, cost is probably one of the main issue.
Other than that, PacBio sequencing also requires more input DNA material (~3-10 ug for human) while Illumina sequencing works with just 100ng. This makes PacBio sequencing more challenging when samples are rare. I'm sure people at PacBio are trying to resolve some of these concerns already though.


text:
Python string editing in Snakemake
I'm running a custom perl script using Snakemake. I use this rule:
rule complexity_20mer_counter:
    input:
       os.path.join(fastq_trimmed_dir, '{sample}_{read}' + fastq_extension.split('.')[0] + '_val_' + '{read}'.strip('R') + '.fq.gz')
    output:
        os.path.join(custom_qc_dir, '{sample}_{read}_20mer_counter.txt')
    shell:
        '''
        #!/bin/bash
        perl scripts/20mer_counter.pl {input} > {output}
        '''

It fails (I'll get to why later) with this error:
Building DAG of jobs...
MissingInputException in line 1 of /rst1/2017-0205_illuminaseq/scratch/swo-406/2019-01-24/rules/20mer_counter.smk:
Missing input files for rule complexity_20mer_counter:
fastq_trimmed/xs0067_P2018SEQE27R01_S1_R1_001_val_R1.fq.gz

... this is because I use the output of trim_galore, this rule looks like this:
rule trim_galore:
    input:
       [os.path.join(fastq_raw_dir, '{sample}_'+ read + fastq_extension)
       for read in reads]
    output:
        timmed_fastq_files = [os.path.join(fastq_trimmed_dir, '{sample}_' + read + fastq_extension.split('.')[0] + '_val_' + read.strip('R') + '.fq.gz')
        for read in reads],
        trimming_reports = [os.path.join(trim_galore_qc_dir, '{sample}_'+ read + fastq_extension + '_trimming_report.txt')
        for read in reads]
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        trim_galore --fastqc --gzip -o {fastq_trimmed_dir} --paired --retain_unpaired {input}
        # Move all qc reports to the fastq_trimmed_dir
        mv {fastq_trimmed_dir}/*.zip \
           {fastq_trimmed_dir}/*.html \
           {fastq_trimmed_dir}/*.txt \
           {trim_galore_qc_dir}
        '''

As you can see, trim_galore not only uses my {read} parameter in the output filename but again adds it at the end, only now without the "R".
I tried solving it in my first rule by applying a python .strip('R') on my {read} wildcard, but this does not seem to work. How can I correctly get '1' and '2' from 'R1' and 'R2' (which are my {read} wildcards)?
(Extra bonus points for a nice output.trim_galore.something as input for my complexity_20mer_counter rule, which does not work because the output is a list (it needs 2 fastq files at the same time), but my complexity_20mer_counter rule can be run in parallel again.)

For me it's quite hard to give you a final solution, without having a look on the whole Snakefile or at least a runnable minimal example. I would suggest that the values for read are already without the R. Than you can easily decide when a R is needed in the filename and when not.

This Snakefile is quite difficult to follow, and its complexity means that subtle mistakes might easily creep in. For example, in your trim_galore rule you have timmed_fastq_files instead of trimmed_fastq_files as the label for your output files. I don't know if this is a problem, but little things like that can definitely mess up the whole pipeline. Also, the shebang at the beginning of each "shell" directive is unnecessary. :-)

Yes, I already caught the timmed error, it was not it. I did try putting reads to ['1', '2'] and making _R part of the string but may people use my pipelines and I want to make the default what every is used most. Doing what @finswimmer suggests will lead to a lot of code to make the "R" optional somehow. Although one could probably figure something out. I would just like to do string operations in the wildcard, like {read}.strip('R'). I'll try to look into making something minimal working but that is a lot of effort in this large pipeline.

It seems that it is not possible to apply strip() directly on the wildcard placeholders. Luckily it is possible to executed functions in the input directive.
So try changing your input of the complexity_20mer_counter rule to this:
lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

The wildcard object is passed to our anonymous function. So have access to the wildcards value via wc.sample and wc.read, which are now handled as normal strings and we can apply strip on them.

Yes, this is what I was looking for, it works, thank you! (I did call "wildcards", not "wc", not sure if that is supported, I call "wildcards" everywhere.


text:
Running Snakemake in one single conda env
I have been experimenting a lot lately with Snakemake, I love it. Recently I also switched to using conda (--use-conda) in the way that is advertised. However, I have some issues with it, mostly related to the way we work. I work in a lab that processes data from our sequencer but also external datasets, moreover, multiple people use my pipeline.
So, with respect to conda this means: currently I have 5 envs, averaging 1.5 GB per env, this takes about 20 min to set up and it's quite some data to pull down. Especially since it is pulled down every time someone analyses data (sometime multiple times a day). This feels like quite a waste of resources. Also, I fix some package versions in the yaml files, but certainly not all (this may lead to issues when repos get reshuffled I read), so I also can't guarantee the exact environment long term (say over years, dependencies may change/update).
So I'm considering switching to a central conda env which I version and of which I export lists of all packages installed for reproducibility reasons. This will then contain all tools.
My question is, is there a nice way to do this? I had issues with putting a "source activate" line in every shell command (and how to do this for Python code under "run" anyway) before, see this: https://stackoverflow.com/questions/54024862/snakemake-conda-env-parameter-is-not-taken-from-config-yaml-file.
So, is there a way, also when using a qsub cluster, to always use a specific conda env, for everything you do?

Using only the --use-conda directive, snakemake copies environments somewhere in the local .snakemake directory. This is repeated for each directory you run it in and replicates all environments effectively for each dataset.
You can however tell snakemake to use a shared environment location using --conda-prefix as described here:

Specify a directory in which the ‘conda’ and ‘conda-archive’
  directories are created. These are used to store conda environments
  and their archives, respectively. If not supplied, the value is set to
  the ‘.snakemake’ directory relative to the invocation directory. If
  supplied, the –use-conda flag must also be set. The value may be given
  as a relative path, which will be extrapolated to the invocation
  directory, or as an absolute path.

This way you could only download/copy the environment once and share it, e.g. this way:
snakemake --use-conda --conda-prefix /path/to/shared/environment

Yes, this is what I was looking for, great, thanx. In this way I can still create multiple environment but they are created once and available to everybody.

--conda should be replaced by --use-conda. Former is not supported by snakemake as far as I could tell.

Your are absolutely right, this has probably changed in new versions. I've replaced it and also the link was updated. Thanks for the pointer

A good practice is to provide yml file along with your Snakefile. This file should contain all the enviroment definitions that are needed to run your pipeline. The minimum content of this file would be to install snakemake itself. Any other program used in your pipeline should be included as well. For reproducibility version numbers must be used.
So, if someone would like to use your pipeline, (s)he first creates a new enviroment based on your definition:
$ conda env create -f my_pipeline.yml

Afterwards activate this environment and run the pipeline:
$ conda activate my_pipeline
$ snakemake

The conda directive in your Snakefile is only neccessary if you need to use program which dependencies cannot be resolved within the main environment. If you start snakemake with the --use-conda parameter, snakemake will create a new environment based on the provided definition file and activate it for the rule.

Especially since it is pulled down every time someone analyses data (sometime multiple times a day).

The environment creation is only done once. If you run the pipeline a second time snakemake can activate the environment diretcly without downloading and installing something again.
The tutorial I wrote some time ago on biostars about snakemake and conda might be useful for you.

Ah, sounds good, I see that "run:" parts also don't support a conda env, so you really need the main env to have everything as well. I'm testing this now, thanx. Will mark as accepted answer when everything works.

Inside your Snakefile you can provide the same conda: directive to each rule. If you install your pipeline in a central location and have everyone use the same conda when running it (hint: use a wrapper script that sets the path to conda) then everyone will be using the same environment. This has the added benefit that you can make this readable to everyone, but disallow them to make any modifications (so users can't mess things up).
The tack we've taken is to have a basic conda env with snakemake and little else and then separate envs for various pipelines (and versions thereof if software versions in your YAML files change) or parts thereof, in cases where you have rules using tools with conflicting dependencies. This ends up being fairly convenient for more complicated pipelines, though if you can get away with just a single environment then finswimmer's answer is simpler.
Note that all of this should work with any scheduler, since snakemake will insert the environment activation stuff in the script it submits to your scheduler.
As an aside, this is exactly how snakePipes works.


text:
Seeking explanation of the hg38 files downloaded from bowtie 2 website
I downloaded the H. sapiens, NCBI GRCh38 files from Bowtie's website. After unzipping, there are 6 files, 4 that end in set 1.ebwt, set2.ebwt, set3.ebwt, and set4.ebwt and two that end in set.rev.1.ebwt and set.rev.2.ebwt.
I am unable to find any explanation for what these files are either online.
which one do I use for bowtie2 mapping?
Thanks

Bowtie is providing you with the index files for hg38.  This are the result of the bowtie2-build indexer.  However you have downloaded the indexes for Bowtie1 not bowtie2

Bowtie 2’s .bt2 index format is different from Bowtie 1’s .ebwt
  format, and they are not compatible with each other.

Then according to the manual you pass the basename of the index for the reference genome. So this should be the name prior to the extension, .1.bt2 for Bowtie2.
bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i>} -S [<sam>]
Main arguments
-x <bt2-idx>
The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final .1.bt2 / .rev.1.bt2 / etc. bowtie2 looks for the specified index first in the current directory, then in the directory specified in the BOWTIE2_INDEXES environment variable.


text:
Building RNAz on MacOS
Having trouble building RNAz 2.1 on MacOS. Following instructions at https://www.tbi.univie.ac.at/software/RNAz/
tar -xzf RNAz-2.1-tar.gz
cd RNAz-2.1
./configure
make

make fails with
/usr/bin/clang++ -o RNAz RNAz.o rnaz_utils.o svm_helper.o zscore.o cmdline.o strand.o svm.o -L../librna -lRNA -lm
Undefined symbols for architecture x86_64:
  "_HairpinE", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
  "_LoopEnergy", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
ld: symbol(s) not found for architecture x86_64

Thinking that compiling for 32-bit might work I tried setting
CXXFLAGS = -m32
LDFLAGS = -m32

without success (same result). Same result with gmake.
Configure
Configure checks for gcc
...
checking for gcc... gcc
checking whether the C compiler works... yes
...

Makefile
Makefile specifies the following:
CC = gcc
CCDEPMODE = depmode=gcc3
CFLAGS = -g -O2
CPP = gcc -E
CPPFLAGS =
CXX = g++
CXXDEPMODE = depmode=gcc3
CXXFLAGS = -g -O2

Compiler
Apple LLVM version 10.0.0 (clang-1000.11.45.5)
Target: x86_64-apple-darwin18.2.0
Thread model: posix
InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin
Has anyone built this successfully on MacOS? Or perhaps suggestions as to how I might build this successfully? I'm at a loss.

Compile with gcc or remove inline from librna/fold.c:776 and librna/fold.c:846.

The other option would be to use RNAz via a docker container instead on installing it directly.  The docker container can be downloaded from biocontainers.  If you are new to docker here is a tutorial for mac.

The first thing to try is to perform the install as sudo.
sudo ./configure # definitely
sudo make ... sudo make install # possibly (?)

The clang++ compiler may be using a path that is not setup in your profile and sudo will overcome this. If sudo works you know what the problem is and its a case of changing permissions and configurations accordingly. There is a good chance this will resolve the matter.
There are other possibilities, it could be the compiler. It looks like the compiler is fine however looking at the config file and the check. The final issue is it could be the program ;-)

Tried this. Same result. Will continue scratching head

Try an installation on Linux
Otherwise its usual stuff, namely upgrading Xtools

Tried this. Same result. Will continue scratching head.

Xtools or Linux ... Linux has a better chance

Xtools is current. Yeah. I know things will likely compile A-OK with Linux. Will try that too. Thanks.


text:
How to model the effect of acid/base on pH
I would like to model a biological compartment with sysBio R package. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:
addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")

which caused continuously decreasing pH levels. I'll leave out the intermediate steps but the simplified version of the last code looks like this:
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "pH", 7.4)
addSpecies(glu, "Acid", 0)

addMAreaction(glu, "Glucose = Lactate + 0.01*Acid", "forwardL",
              "backwardL", name = "glycolysis, sort of")

addMAreactRate(glu, "forwardL", "assigned", "kgl*Glucose")
addMAreactRate(glu, "backwardL", "assigned", "ka*pH+klg*Lactate")

addParameters(glu, "kgl", 0.01)
addParameters(glu, "klg", 0.01)
addParameters(glu, "ka", 0.05)
addParameters(glu, "a", 0.001)

addRule(glu, "pH rule", "ODEs", "pH=7.4-a*Acid")

makeModel(glu)
simResults <- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")

And the resulting plot is this where pH is increasing forever despite Acid is at a steady state concentration.

Please note that the abovementioned code is just a simplification. It doesn't even look like aerobic glycolysis where 1 glucose produces  lactate and 2 hydrogen atoms etc.
I intend to include several acid- and base-creating reactions and evaluate the resulting pH. So could someone point me in the right direction about how to model the relationship between pH and acids or bases?
The answer does not need to be explicit code. I'll appreciate any solid comment on my design errors.

I'm unfamiliar with sysBio, but it looks like you may just be plotting a relationship, and need to know the relationship. Is the Hendersen Hasselbalch equation running behind the scenes, or is that what you're looking for? I'd note, if you're trying to model acid base in human physiology, you're making a number of assumptions by leaving out the major systems for acid base regulation.

You are right. But please keep in mind that this is just a minimal working code sufficient to reproduce the error. You may think of it like a reaction where A reversibly converts to 2*B + 2+H, and H reversibly causes P to sink to 0. At least this was how I modeled it. For example I tried to add rules to limit pH within 0 and 14 as well but these are not included in the example. Apparently some familiarity with either SBML, ordinary differential equations, or maybe sysBio is needed. I couldn't reach the original author of the package so I'm trying to see if I'm missing sth fundamental.

Since your question pertains to the use of a package and not really about the biology behind the model, your question would not be on-topic here. It should, however, be on-topic on [bioinformatics.se].

I figured that one has to include only the changes in the rules formula. Keeping the initial value of the Species or Parameter when defining circular or self-referencing formulas caused the error. So the correct rule becomes:
pH=-a*Acid instead of pH=pH-a*Acid or pH=7.4-a*Acid.
For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL",
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP",
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO",
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned",
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned",
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned",
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned",
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned",
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned",
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

The head of the resulting data.frame is:
  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

and resulting graph looks like:


text:
Hardware Requirements (specs) for Bioinformatics-dedicated desktop
I know this is a somewhat general or vague question, but I’m interested in your opinions. I must build a desktop for general bioinformatics activities with human genomes. I will work with Python and R libraries, I will also need to install denovo sequence assembly software, alignment tools, variant viewers/browsers, sql, among other software. I don’t have much experience customizing hardware that will be dedicated to bioinformatics and genome data analyses. I’m pretty sure I work with linux and/or iOS. What is the recommended RAM, processor speed, HD capacity, graphics cards, etc. I appreciate your opinions on the matter. Thank you.

You mention sequence assembly software. For the specs of the machine it is important to know what you will be assembling. Are you assembling genomes or transcriptomes? And will you be assembling human genomes/transcriptomes or other organisms?

The RAM requirements for some bioinformatics analyses like assembly can be quite high (in the hundreds of Gigabytes). My recommendation is to get a fast laptop. Something with an i7 quad-core processor, 16 GB of RAM, and 1TB of storage should do. Do your analyses on this laptop. If you run out of RAM or storage, then use a cloud service to get access to bigger machines or buy an external hard drive. Once you have a feeling for what your actual needs are, start looking at a dedicated rack server solution in your local institution's server farm or again, in the cloud.
You don't mention any analyses that clearly benefit from a GPU, but as you get more familiar with your analyses, perhaps that will become useful as well.

Thank you. I believe your advice is quite reasonable, especially since I’m not sure yet what species I’ll be working with. In all likelihood, I’ll have access to server based software for large sequences, but the more I can do locally the better.

I am glad you feel that way. If you think this answer is a good one and is the best answer to your question, then I recommend you accept this answer.

Yes. I realize how vague my question is. I don’t have much information on what resources I’ll have access to yet, but I was asked to give some specs so my PC can be configured accordingly. My best answer at this point is to be provided with something that can easily be upgraded if needed. Thanks again.

If you deal with genomics NGS data then heavy RAM will be essential. I don't think you will get this from a reasonable priced laptop.
If you are using a Deep Learning based algorhithms or Bayesian or maximum likelihood heavy duty processors are needed.
I would go for an allround machine.
Two hard drives:

0.5 to 1 Tb SSD (I use 1TB),
2 to 4 Tb Stata data storage,
4-core (light) or 8-core, i7 (10-core is common now),
16GB-32GB RAM.

Important: For your purpose if you are doing de novo NGS assemblies use:

ALOT of RAM.
Sufficient harddrive - make sure this is backed up

If you want evidence why RAM is important simply look at the threads here and on Stackexchange which read ... de novo NGS assembly PC has stopped working?(!).
The following threads give an insight into the sort of problems RAM creates

here:
a similar RAM problem is here, but not related to NGS is here

The following looking like a particularly bad case:

grasshopper NGS


text:
How can I interpret gene expression data from Bioconductor packages?
I am currently looking at microarray data from a bioconductor microarray dataset. Specifically, I have data (a snippet) which looks like the following:
   GSM794885   GSM794886   GSM794887   GSM794888   GSM794889
   -1.67170290 -0.07803822  0.99990560 -1.94692850 -1.78671360

I am wondering how to interpret the values in this microarray dataset. Is it generally the case that larger values mean a more significant effect?

Try to read more about microarray data analysis first, follow some of the tutorials, and you will understand what these values are.

Welcome to the site. Please follow benn advice! Also this values should be significant for what? Provide more information about how did you obtained these values (what is the code needed to reproduce this values)?

Looks like a log2(intensity value), but without more information it's impossible to say for sure...

It would be great if you could add more context to your question for those who want to help you out with your problem. What type of microarray is this? What are you trying to do with this dataset? What question are you trying to answer?

Processed microarray data is commonly represented as log_2(expression). This data transformation is used because the data more closely fit a normal distribution in log space. With such a transformation, negative values represent very small amounts, rather than the absence of something.
Larger values mean that more stuff is being produced (presumably transcripts), but its significance can only be assessed when combined with other data (e.g. through differential expression analysis).


text:
Perfect Phylogeny vs Maximum parsimony
I am searching various sources about phylogenetics. I saw some materials about perfect phylogeny and also phylogenies acquired from maximum parsimony constraint. They seem very similar to me. Are they the same?

No. A perfect phylogeny is one such that all characters evolve on the tree with no homoplasy, i.e. for a binary character, changes occur once from 0 -> 1 but never from 1 -> 0. Maximum parsimony inference may produce a perfect phylogeny, but typically for real datasets some degree of homoplasy is required to explain character patterns.

Long discussion that took many years to resolve.
MP doesn't account for back mutation, which is a HUGE problem for nucleotide data because in theory 1:4 mutations is a back mutation. Maximum likelihood has resumed is crown here, however Beast deployment of a Bayesian calculation is also very widely, particularly for molecular dating. ML uses  a reversible matrix A <-> T, whilst Beast uses a directional matrix A->T and will score T->A separately.
Happy to discuss the approaches and modern interpretations at length.
The model which accounts for back-mutation is the Jukes-Cantor correction (JC correction), this is a basic model present in every phylogenetic algorithm except p-distances). Essentially JC extrapolates back-calculation for the true number of mutations against the observed number of mutations against divergence time. However, when the P-distance (uncorrected observed distance) exceeds 0.75 for nucleotides the model is not viable. Basically, nucleotide divergence >0.75 is saturated and the phylogenetic information is essentially random.

Imo a general comment about phylogenetic methods is not particularly helpful for the distinction between most MP trees and one that mathematically fits the description of a perfect phylogeny. But to which models are you referring to for ML and BEAST? AFAIK they have consistent sets of CTMC models.


text:
Error in installing Vep
I am trying to installing Vep on Linux but an error is permanent
[fi1d18@cyan01 ensembl-vep-release-95]$ perl INSTALL.pl

Hello! This installer is configured to install v95 of the Ensembl API for use by the VEP.
It will not affect any existing installations of the Ensembl API that you may have.

It will also download and install cache files from Ensembl's FTP server.

Checking for installed versions of the Ensembl API...done

Setting up directories
Destination directory ./Bio already exists.
Do you want to overwrite it (if updating VEP this is probably OK) (y/n)? y
 - fetching BioPerl
 - unpacking ./Bio/tmp/release-1-6-924.zip
 - moving files
Attempting to install Bio::DB::HTS and htslib.

>>> If this fails, try re-running with --NO_HTSLIB

 - checking out HTSLib
fatal: destination path 'htslib' already exists and is not an empty directory.
 - building HTSLIB in ./htslib
In /home/fi1d18/ensembl-vep/ensembl-vep-release-95/htslib
gcc -g -Wall -O2 -fPIC -Wno-unused -Wno-unused-result -I.  -c -o cram/cram_io.o cram/cram_io.c
cram/cram_io.c:61:18: error: lzma.h: No such file or directory
cram/cram_io.c: In function 'lzma_mem_deflate':
cram/cram_io.c:679: warning: implicit declaration of function 'lzma_stream_buffer_bound'
cram/cram_io.c:685: error: 'LZMA_OK' undeclared (first use in this function)
cram/cram_io.c:685: error: (Each undeclared identifier is reported only once
cram/cram_io.c:685: error: for each function it appears in.)
cram/cram_io.c:685: warning: implicit declaration of function 'lzma_easy_buffer_encode'
cram/cram_io.c:685: error: 'LZMA_CHECK_CRC32' undeclared (first use in this function)
cram/cram_io.c: In function 'lzma_mem_inflate':
cram/cram_io.c:695: error: 'lzma_stream' undeclared (first use in this function)
cram/cram_io.c:695: error: expected ';' before 'strm'
cram/cram_io.c:701: error: 'LZMA_OK' undeclared (first use in this function)
cram/cram_io.c:701: warning: implicit declaration of function 'lzma_stream_decoder'
cram/cram_io.c:701: error: 'strm' undeclared (first use in this function)
cram/cram_io.c:701: warning: implicit declaration of function 'lzma_easy_decoder_memusage'
cram/cram_io.c:719: warning: implicit declaration of function 'lzma_code'
cram/cram_io.c:719: error: 'LZMA_RUN' undeclared (first use in this function)
cram/cram_io.c:720: error: 'LZMA_STREAM_END' undeclared (first use in this function)
cram/cram_io.c:732: error: 'LZMA_FINISH' undeclared (first use in this function)
cram/cram_io.c:743: warning: implicit declaration of function 'lzma_end'
At top level:
cc1: warning: unrecognized command line option "-Wno-unused-result"
make: *** [cram/cram_io.o] Error 1
Compile didn't complete. No libhts.a library file found at INSTALL.pl line 820.
[fi1d18@cyan01 ensembl-vep-release-95]$

And
[fi1d18@cyan01 ensembl-vep]$ perl INSTALL.pl --NO_HTSLIB

Hello! This installer is configured to install v95 of the Ensembl API for use by the VEP.
It will not affect any existing installations of the Ensembl API that you may have.

It will also download and install cache files from Ensembl's FTP server.

Checking for installed versions of the Ensembl API...done

Setting up directories
Destination directory ./Bio already exists.
Do you want to overwrite it (if updating VEP this is probably OK) (y/n)? y
 - fetching BioPerl
 - unpacking ./Bio/tmp/release-1-6-924.zip
 - moving files

Downloading required Ensembl API files
 - fetching ensembl
 - unpacking ./Bio/tmp/ensembl.zip
 - moving files
 - getting version information
 - fetching ensembl-variation
 - unpacking ./Bio/tmp/ensembl-variation.zip
 - moving files
 - getting version information
 - fetching ensembl-funcgen
 - unpacking ./Bio/tmp/ensembl-funcgen.zip
 - moving files
 - getting version information
 - fetching ensembl-io
 - unpacking ./Bio/tmp/ensembl-io.zip
 - moving files
 - getting version information

Testing VEP installation
ERROR: Testing VEP script failed with the following error
Can't locate Time/HiRes.pm in @INC (@INC contains: /home/fi1d18/ensembl-vep/modules /home/fi1d18/ensembl-vep ./Bio /usr/local/lib64/perl5 /usr/local/share/perl5 /usr/lib64/perl5/vendor_perl /usr/share/perl5/vendor_perl /usr/lib64/perl5 /usr/share/perl5 .) at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Utils/SqlHelper.pm line 93.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Utils/SqlHelper.pm line 93.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBConnection.pm line 81.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBConnection.pm line 81.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBAdaptor.pm line 69.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/DBSQL/DBAdaptor.pm line 69.
Compilation failed in require at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Registry.pm line 137.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/Bio/EnsEMBL/Registry.pm line 137.
Compilation failed in require at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseVEP.pm line 58.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseVEP.pm line 58.
Compilation failed in require at (eval 7) line 3.
        ...propagated at /usr/share/perl5/base.pm line 94.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/BaseRunner.pm line 56.
Compilation failed in require at (eval 6) line 3.
        ...propagated at /usr/share/perl5/base.pm line 94.
BEGIN failed--compilation aborted at /home/fi1d18/ensembl-vep/modules/Bio/EnsEMBL/VEP/Runner.pm line 71.
Compilation failed in require at ./vep line 20.
BEGIN failed--compilation aborted at ./vep line 20.

[fi1d18@cyan01 ensembl-vep]$

Do the VEP installation instructions give a list of dependencies or prerequisites? If so, you must install these before trying to install VEP. The solution to your second error is also in the message: Can't locate Time/HiRes.pm. This means there's a Perl module Time::HiRes that is not installed.

A general tip with Perl installations: When you see an error message that states "Cannot locate <thing1>/<thing2> in @INC", the solution is almost always running cpan <thing1>::<thing2> from the command line.

Says You need to be root to perform this command

It is pretty rare for linux utilities to need root. Try Google "CPAN install module no root" or "CPAN user install module". By adding a "without root" to your Google searches, you will see ways to not need root.

Most "install" processes copy files and compile them at a certain location (ideally something found under $PATH), creating executable files so when you run the command, the shell can find the executable file. When you don't have root, you might be unable to copy the files to the default location, but you would be able to do the same set of install steps to a custom location, which you can then add to $PATH yourself. Once you understand the concept of how installs work, you'll seldom need root access.

What you can try is stated out in the message:

If this fails, try re-running with --NO_HTSLIB

Thank you, I got error again. Please have a look at my error when I am trying with --NO_HTSLIB


text:
wget for links inside html pages
I am trying to download a file from the following repository: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474
As you can see, there are several layers to the webpage. For example, clicking on the download tab doesn't change the URL and, the link is not a 'download link' per ce - where simply clicking the link automates a download. I have tried some of the answers on the forum, where they have advised using quotes and the operation:
wget "url/?target=link". This however does not work in the following instance.
Basically, I want to download the file labeled 'P1TLH.bam', in the download tab of the link provided:

Thanks, found it. But doesn't wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam work?

It does! Thank you so much. What is the logic behind it though? Don't quite understand why the original command was not working..

Because you weren't using the link to that file, presumably. Wget will download what you tell it to, so you need to tell it to get the target of the link you want.

If you want to download a file, you need to use the link to that file. Your original attempt, wget  https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474 wouldn't work since that's a link to the trace page of the relevant run. If you want to download something else, just right click on the link (the one in the screenshot in your question), copy the URL and use that:
wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam

NCBI has its own digger for this purpose. No, my bad the traditional Perl approach of LWP and HTTP::Request::Common to interrogate NCBI, isn't going to work no how much data munging or spider crawling is done, because its just an NGS data dump. ... well its diggable, but it would be easier to dump the data and  interrogate it locally.

Please post the link and an explanation of how to use it, yes. Without it, this isn't really an answer. With it, it would become a very useful answer, much more so than my own.

Well it obviously isn't! :0) Its in Perl, so no-one understands that now

Oy! Nothing wrong with Perl! Just because these young whippersnappers don't get its beauty, doesn't mean it isn't a great language! Joking apart, please edit the answer so it can stand on its own or just delete and post a comment instead. As you see, it's already attracted a downvote (not mine).

No, I just didn't get around to doing it and when I get around to it I'll update the post

authors: محمدرضا افشار


text:
Parallel version of cn.mops' bamDataRanges
I am trying to run cn.MOPS to detect CNVs in multiple bam files. Following the example code, I have run:
source("https://bioconductor.org/biocLite.R")
biocLite("cn.mops")
options(width=75)
set.seed(0)
library(cn.mops)
library(Biobase)
library(GenomicRanges)
library(GenomeInfoDb)
cn.mopsVersion <- packageDescription("cn.mops")$Version
library(cn.mops)
BAMFiles <- c('bam1.bam', 'bam2.bam', ..., 'bamN.bam')
bamDataRanges <- getReadCountsFromBAM(BAMFiles)

PLEASE BE PATIENT... this might take a while. Consider using the parallel version of this function

Since this is indeed taking a while, I'd like to follow the tool's suggestion and use "the parallel version of this function". But how do I do that? I don't see any mention of "parallel" in the manual. Is there some standard Bioconductor/R trick to make a function parallel? Am I missing something in the manual?
How can I run bamDataRanges parallelized?

It looks like parallel is an argument you set and feed to the function? https://rdrr.io/github/Bioconductor-mirror/cn.mops/man/getReadCountsFromBAM.html

@Chris_Rands duh! Thanks, I should have found that myself. Could you post it as an answer?

@Chris_Rands I posted an answer just so the question can be counted as answered. Please feel very free to post your own (even just copying mine) and I'll delete and accept yours.

This is an option for the function documented in its help page. Just run ?bamDataRanges and you will see:
Usage:

     getReadCountsFromBAM(BAMFiles, sampleNames, refSeqNames, WL = 25000,
       parallel = 0, ...)

[...]

parallel: The number of parallel processes to be used for this
          function. Default=0.

So all you need to run with 10 cores is:
bamDataRanges <- getReadCountsFromBAM(BAMFiles, parallel=10)


text:
Calculate the percentage of each unique phylogenetic tree in a BEAST output
I have a nexus formatted BEAST output containing 20,000 phylogenetic trees of seven taxa. Is there any way to get the percentage of each unique phylogenetic tree contained in this output?
I already made an unsuccessful attempt with R.

If the phylogenetic trees can be represented as strings then it is basically just string matching. Anyway, I think your question would be more suited for [bioinformatics.se] as it is mostly about the programming aspect.

Could you explain what is a "unique phylogenetic tree"? Perhaps if you upload an image it would be clearer what do you mean.

@llrs When I mention a 'unique phylogenetic tree', I am referring to the topology of the tree. I would like to extract all the unique topologies from my BEAST output, and calculate how many times each of these topologies appears in my BEAST output. This would allow me to establish a percentage for each topology. The idea is to know which topologies are the most represented in addition to the consensus tree.

I finally managed to do it in R. Here is my code:
install.packages('devtools')
library(devtools)
install_github('santiagosnchez/rBt')
library(rBt)

beast_output <- read.annot.beast('beast_output.trees')
beast_output_rooted <- root.multiPhylo(beast_output, c('taxon_A', 'taxon_B'))
unique_topologies <- unique.multiPhylo(beast_output_rooted)

count <- function(item, list) {
  total = 0
  for (i in 1:length(list)) {
    if (all.equal.phylo(item, list[[i]], use.edge.length = FALSE)) {
      total = total + 1
    }
  }
  return(total)
}

result <- data.frame(unique_topology = rep(0, length(unique_topologies)),
                     count = rep(0, length(unique_topologies)))
for (i in 1:length(unique_topologies)) {
  result[i, ] <- c(i, count(unique_topologies[[i]], beast_output_rooted))
}

result$percentage <- ((result$count/length(beast_output_rooted))*100)

This is fantastics, thanks for posting back

You require a consensus phylogeny.
This is available in RAxML, PAUP, the old Phylip package and possibly the stats companion to Beast, Tracer. I'd need to dig up the code for RAxML and PAUP. The easiest is Consense here because you just download it and input the sequence. NOTE their is possibly a format difference (RAxML and Consense use Phylip format), I suspect you've got nexus, or Newick format.
Beast is likely to have its own consensus formating OR the Beast companion  "Tracer" may have this (I can't remember). Tracer is here and would be cool if it did consensus phylogeny because there will be no conversion problems (it doesn't use Phylip format). The alternative would be to use PAUP (no conversion problems), which is currently in a state of flux from various stages of freeware and updates. PAUP is contained within the Genious package, which is available as a 1 month free-ware.
Have a look and let me know: because if you want to use RAxML or Consense you'll likely need to remind me of the Beast output, if it is what I think it is a format conversion app is needed.  The calculation is easy, its easy to script a conversion program, but I think tree formating will be an issue.
You could use Ape, within R, but as you have gathered there are alot of other solutions. There are more recent, impressive, solutions as stand-alones which I've not used yet.

Many thanks for your answer! However, I already have a consensus tree and I am now trying to get the other most represented topologies from my BEAST output.

Ok. As far as I'm aware that program has not been written, at least a couple of years ago it wasn't, not sure about this now. Using R "ape" is your best bet. On this question I definitely can't help.


text:
General framework for fusing biological networks
I know that there are different types of biological networks. Protein-protein interaction, gene interaction, metabolic networks, gene expression networks, etc. Each of these networks can be represented as a graph.
1) Is there a principled way to integrate these graphs into a single large graph? For example, probably most of these networks have nodes that are genes, or gene related things, and this could allow us to match nodes in one network to nodes in another network.
2) Is there a principled way to ask joint questions about these networks?

Hi Mendel, welcome to Bioinformatics StackExchange. If you want more informative answers, it would be great if you could add some more context around your question (i.e. what particular problem are you trying to solve; be as specific as possible). There are a lot of different ways to approach bioinformatics problems, and this leads to a lot of garden-path discussion that doesn't really end up helping the person who has asked the question.

I believe you're referring to something related to multi-layer networks. Here you can find an interesting reading about this kind of networks:
https://academic.oup.com/comnet/article/2/3/203/2841130
Also you can try with this one:
https://www.sciencedirect.com/science/article/pii/S0370157314002105
As a matter of fact, I'm not an expert of multi-layer networks but, as far as I can tell, the main problem (as you mentioned) is to find correlations, or connections, between the different layers in order to be able to draw some meaningful conclusion.
I am not sure you can call this a multi-layer network but this can be an example of what you're looking for: suppose you have a protein-protein interaction network. You can enrich your network using biological information. You can associate biological functions, or cellular locations or whatever you can find in GeneOntology (just to make an example of a biological database), as a further layer of information regarding your protein-protein interaction network.
This will help you in elucidating the role a node (or more nodes) is playing with respect to its neighbours: let's say five nodes form a clique and you also know they share the same cellular location...then you're able to say something about these nodes (maybe they form a complex?) that has a dual nature: biological and graph-based.
Hope it helps!

You are talking network statistics and this is a complex area of stats.
For generalised simplicity, Google fusion tables have a simple to use but impressively interactive network tool. This is a great starting point for exploration. Have a go and see what you think.
If you are looking at a haplotype map, there is a formal solution to this (parsimony network), which is standalone program written some time ago but still useful today. In my area split-decomposition is an important important technique.
The easiest question is the number of connecting nodes and Google fusion tables is as good as any routine interactive networking tool. As I say beyond this network stats is a specialised field.


text:
Save output to a specific folder and/or with a specific prefix in Cancer Genomics Cloud
I would like to run the CNVkit workflow on Cancer Genomics Cloud.
However, the two workflows that I'd like to use:

CNVkit Reference, for building a reference from normal samples
CNVkit CNV Calling, for calling CNVs by comparing tumors to the reference

don't have parameters that can be set to define in which folder to save the results or which prefix to add to the results. This means that if you're generating multiple references or calling CNV for multiple tumor samples you will end up with results having the same name for each run (tumor).
I tried to modify the CWL workflow, online and through Rabix composer, but I was unable to make it work. I would like to use the published workflow to not code everything from the ground up.
Do you know how to update/use the current CWL workflow in order to avoid having the same filename/folder for all the results?

@llrs updated question with links to the workflows' code

If the output is File, you can add this step https://github.com/Barski-lab/workflows/blob/master/tools/rename.cwl to rename the output file.


text:
Percentage distribution of cells in all clusters based on their treatment condition?
I have 2151 cells, I clustered them by Seurat to 5 clusters. With the code below, I am able to have the number of cells per cluster and per condition:
number_perCluster<- table(object@meta.data$conditions,
                          object@meta.data$clusterID)

Could I please get a hint on how I could proceed to produce a figure where I say: X% of cells treated with condition "Y" are located in cluster "C" and so forth?
Thank you in advance.

Have a look at prop.table, also if you want a figure, please describe it a bit more. This is usually better displayed as a table or heatmap, but perhaps you want something else...

Please don't also create a github issue on the Seurat page each time you ask a question here, there is no need

Here is a solution using dplyr and ggplot2:
library(Seurat)
library(dplyr)
library(ggplot2)

meta.data <- pbmc_small[[]]

# create random classifications for the sake of this example
meta.data$condition <- sample(c('A', 'B', 'C'), nrow(meta.data), replace = TRUE)

counts <- group_by(meta.data, condition, res.1) %>% summarise(count = n())

ggplot(counts, aes(res.1, count, fill = condition)) +
  geom_bar(stat = 'identity')


text:
sva for RNA-Seq data without known phenotype
I have been working on RNA-Seq data from two different cohorts, and they show very strong batch effect (~35% variance explained by 1st component in PCA). Since I am trying to do a class discovery from a data set with the subtype of only some samples are known, the only methods I have been using are ComBat and pSVA from SVA package. I would like to compare them with more methods to make the clustering result more trust-worthy.
I found that some papers used sva package, but not citing the paper on pSVA, meaning perhaps there are some other ways to use sva for this type of data that I am not aware of. Is there a method to construct a null model and a full model for sva only by knowing the batch the samples are from?
Any other suggestion of batch effect correction tools for a data set with no known phenotype/subtype are also welcome.

You can just remove the first component from the dataset by setting the first eigenvalue to 0 in the diagonal matrix and then multiplying the SVD matrices.
I am not sure which R code you are using to estimate the components (I have a comparison of PCA code in R on https://github.com/aedin/ODSC_2018/tree/master/PCA_vignette)
In ade4 R package, there is a function reconst() which reconstructs data from a selection of components: https://rdrr.io/rforge/ade4/man/reconst.html
Happy to provide code to do this in other PCA versions if needed
Aedin Culhane

Thank you for your suggestion. I am more of a molecular biologist so I don't quite know how to do it. I understand vaguely the concept but not the actually steps. Or is it what ade4 does to remove the first component?