# Question
Title: Calculating celltype proportion changes between conditions in single-cell data
I want to test for significance between differences in cell count between conditions.

I have 2 conditions (N, T), 3 donors (samples) each, and four cell types (A-D) as shown in the example below.

I believe this problem is analogous to the differential expression test performed with edgeR or DEseq.

Example data:

```
Count   Sample  Celltype    Condition
555 P1  D   T
3829    P1  C   T
945 P1  A   T
222 P1  B   T
795 P2  D   T
5464    P2  C   T
839 P2  A   T
438 P2  B   T
629 P3  D   T
6426    P3  C   T
647 P3  A   T
457 P3  B   T
8   P1  D   N
5122    P1  C   N
3204    P1  A   N
1   P1  B   N
20  P2  D   N
6111    P2  C   N
2524    P2  A   N
26  P2  B   N
66  P3  D   N
7490    P3  C   N
5272    P3  A   N
3   P3  B   N

```

# Answer
This publication used a negative binomial regression model of the counts of cells at each time-point.

From the methods section: *"For each cell type, we model the number of lineage-labelled cells detected in each analysed mouse as a random count variable using a negative binomial distribution. The frequency of detection is modelled by using the natural log of the total number of cells of that type profiled in a given mouse as an offset. The time point of each mouse (0, 30 or 60 days post tamoxifen) is provided as a covariate. The negative binomial model was fit using the R command ‘glm.nb’ from the ‘MASS’ package. The P value for the significance of the change in labelled fraction size between time-points was assessed using a likelihood-ratio test, computing using the R function ‘anova’."*

The method is available here as an R script.

The input `x` is a dataframe with the following columns:

```
x$variable # celltype
x$value # cell count
x$batch # condition+donor
x$timepoint # condition (e.g. adult, fetal)

```

Montoro, D.T., Haber, A.L., Biton, M. et al. A revised airway epithelial hierarchy includes CFTR-expressing ionocytes. Nature 560, 319–324 (2018) doi:10.1038/s41586-018-0393-7
> 1 votes

---
Tags: scrnaseq, single-cell
---


# Question
Title: WGCNA correcting p value for varying sample numbers in each of the trait columns
I am correcting a data frame of sample comparison using p value. However in some columns I have different sample numbers. Within the data frame I have 4 comparisons e.g B\_R, R\_C, B\_C, C\_C. They each have different sample numbers in the comparisons. I need to create a loop or change the code so I can vary the nSamples numbers for each of the comparisons. I have tried to do this for each comparison separately but to no avail. Does anyone have any suggestions? Thank you. I have tried..

nGenes = ncol(cells)
nSamples = nrow(cells)

MEs0 = moduleEigengenes(cells, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")  ``````````

I have tried to code for the varying Sample numbers using:

MEs0 = moduleEigengenes(cell, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor= cor(MEs, NR_status, use="p")
moduleTraitPvalue$B_R = corPvalueStudent(moduleTraitCor, 32)
moduleTraitPvalue$R_C = corPvalueStudent(moduleTraitCor, 31)
moduleTraitPvalue$B_C = corPvalueStudent(moduleTraitCor, 47)
moduleTraitPvalue$C_C = corPvalueStudent(moduleTraitCor, 31)

sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue$NR_R,
                      moduleTraitPvalue$R_C,
                          moduleTraitPvalue$C_C,
                      moduleTraitPvalue$NR_C, 1), ")", sep = "")
</code></pre>

# Answer
You can use the function `corAndPvalue`. Simply replace the lines

```
moduleTraitCor= cor(MEs, status, use="p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

```

by something like

```
cp = corAndPvalue(MEs, status)
moduleTraitCor = cp$cor
moduleTraitPvalue = cp$p

```

`corAndPvalue` takes automatically into account the number of valid samples for each element of the resulting correlation matrix.
> 0 votes

---
Tags: r, wgcna
---


# Question
Title: How to show bootstrap values on a phylogenetic tree constructed with RAxML
I apologize for this very basic question. With @Michael G.'s help, I have been learning how to use RAxML to build a phylogenetic tree (How to define an outgroup to build a robust amino acid tree) and after running my analysis I obtained several results files. I have both a *RAxML\_bestTree.result* file and a *RAxML\_bootstrap.result* file. How do I "combine" these files so that my tree shows the bootstrap values?

Any help will be greatly appreciated.

# Answer
> 3 votes
Apologies for the delay. The calculation is a majority consensus phylogeny of a maximum likelihood bootstrap, which is superimposed into a maximum likelihood phylogram. That's the tech. speak over ... the command is,

```
raxmlHPC -f b -m PROTGAMMAILG -n output_bootstrap.tre -t RAxML_bestTree* -z RAxML_bootstrap.result

```
* RAxML\_bestTree\* = your maximum likelihood tree without bootstraps (file is probably 1 line long)
* RAxML\_bootstrap.result = your complete list of bootstraps (file is 500 lines long \[or the number of bootstraps you set\])
* output\_bootstap.tre = final result tree, i.e. majority consensus. This file is one line long but will have with all bootstraps attached to each branch.

The calculation will take less than one minute. The code is C and very slick. Note the "RAxML\_bestTree\*" is probably better to write the file name in full, not truncated via the glob (\*). BTW you must specify the model.

**Visualising the tree** You use a Mac ... so the cool tree program is figtree . Make sure when you load it you click the branch option and set it to see your bootstraps, otherwise it will not visualise them by default. You might need to re-root the tree, from memory your command specified an outgroup. It might not be the right one though and figtree will easily re-root it for you. Bootstrap values above 75% are considered robust. If you get stuck, can chat but I'm not around until this time next week (bereavement leave and loads of meetings).

# Answer
> 2 votes
Depending on the search parameter settings RAxML produces several output files among which you normally find:

> **RAxML\_bipartitions.result**: If you used the -f b option, this file will contain the input tree with confidence values from 0 to 100 drawn on its nodes! It is also printed when -f a x have been specified, at the end of the analysis the program will draw the BS support values on the best tree found during the ML search.
>
> **RAxML\_bipartitionsBranchLabels.result**: Contains the same information as the file above, but support values are correctly displayed as Newick branch labels and not node labels! Support values always refer to branches/splits of trees and never to nodes of the tree. Note that, some tree viewers have problems displaying branch labels, they are however part of the standard Newick format.

from https://cme.h-its.org/exelixis/resource/download/NewManual.pdf

If you don't find the files in your output directory you did not specify the above parameters. If a rerun is not feasible you could you use the multinewick file RAxML\_bootstrap.result with raxml or as input for the program splitstree to generate a consenus tree with different rules (majority rule, strict, etc.). However, these options don't guarantee to generate the topology of the maximum likelihood tree.

---
Tags: phylogenetics, phylogeny, software-usage
---


# Question
Title: Where can I get ensembl Medaka genome for RNAseq
I am trying to map a RNA-seq dataset using ensembl genome for medaka fish.

From here http://uswest.ensembl.org/Oryzias\_latipes/Info/Index

when I click on a `Download DNA sequence (FASTA)` I am forwarded here ftp://ftp.ensembl.org/pub/release-98/fasta/oryzias\_latipes/dna/

where I have a huge number of `.fa.gz` files and I am just looking for a single file, that I can use.

While when I am downloading medaka genome from NCBI (and using this technique) http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download I have my ncbi genome in one file.

How can I get an ensembl genome in one file?

# Answer
> 3 votes
Download the ones without a chromosome number in the filename. Those files are the whole genome.

---
Tags: rna-seq, ensembl
---


title: bedtool intersect multiple bed files

text:
# Question
Title: bedtool intersect multiple bed files
I have 7 biological replicate of a normal sample after peak-calling I converted it into bed files and now I want to find what are the common region between those .

So i tried with bedtool intersect

```
my command is this
bedtools intersect -a SRR2920506.filtered.bed -b SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

```

it says " the intersect tool can detect overlaps between a single -a file and multiple -b files (instead of just one previously). One simply provides multiple -b files on the command line."

My doubt is the file that goes to "a" argument is not query which i have to detect overlaps against "b" files as these are all biological replicate. Am i doing the right way?

Secondly I tried with `bedops` but the number of overlaps which is reported is really low when i compare the result which i get from bedtools intersect function.

Any suggestion or help regarding how to intersect multiple would be really appreciated

# Answer
If you want to find the intersection between ALL your bedfiles, you can try `multiIntersectBed` (available since bedtools 2.14.3).

It should work like this :

```
multiIntersectBed -i SRR2920506.filtered.bed SRR2920531.filtered.bed SRR2920466.filtered.bed SRR2920478.filtered.bed SRR2920507.filtered.bed SRR2920532.filtered.bed  > answer.bed

```

To sort your bedfiles by chromosome then by start :

```
sort -k 1,1 -k2,2n myfile.bed

```
> 5 votes

---
Tags: bedtools
---


# Question
Title: Is BWT based aligner suitable for any types of alignment task?
Burrows wheeler transformation based aligner like BWA or bowtie seems a standard alignment tool used many area. I was just wondering if there is a kind of alignment task in which BWT algorithm is not suitable.

# Answer
> 3 votes
BWT based algorithm is just a component of a full-pledge aligner. To that end, it is applicable to most alignment tasks. However, for long-read alignment and metagenomes, we don't often use BWT based algorithms as other methods are faster.

# Answer
> 1 votes
Not sure about the algorithm itself but Heng Li wrote an article about the future of BWA when releasing minimap2. It gives helpful insights about the topic : Heng Li's blog

# Answer
> 0 votes
Surely you know that neither bwa nor bowtie by itself is splice-aware?

---
Tags: sequence-alignment, bwa, bowtie2
---


title: Find indels between two short sequences

text:
# Question
Title: Find indels between two short sequences
I have two sequences, say AAAGCTCGAGG and AAAGCGAGG. I need a convenient tool which shows me insertions or deletions between these, i.e. in this case something like

```
AAAGCTCGAGG
|||||  ||||
AAAGC  GAGG
2bp deletion starting at base pair 5

```

and also insertions such as

```
AAAGCTC GAGG
||||||| ||||
AAAGCTCTGAGG
1bp insertion at base pair 7

```

Does something like that exist? Blastn finds the alignment of course, but only outputs the *number* of gaps, not their *position*. I guess I could parse a pileup, but do not want to create a bam file from one of the sequences first, so samtools seems to be out. Any further ideas?

# Answer
> 1 votes
Smith-Waterman is guaranteed to give you the best alignment given your choice of substitution matrix. For a few such tiny sequences, it's what I would use. Many implementations exist.

# Answer
> 0 votes
Perhaps you could use something like `SearchIO` to parse the output of BLAST, to get the desired alignment attributes: https://biopython.org/DIST/docs/api/Bio.SearchIO.BlastIO-module.html

# Answer
> 0 votes
blastn should be able to output this format of alignment result, using the `--outfmt 0` option

# Answer
> 0 votes
Is there any particular reason why you don't want a BAM output format? If the sequences are sufficiently long, then the CIGAR field from the SAM output of minimap2 could be used:

```
minimap2 -a seq1.fa seq2.fa | samtools view | awk '{print $6}'

```

---
Tags: samtools, indel, blastn
---


# Question
Title: how to search NCBI to find all available reviewed accession numbers for a protein family
what is the most optimal keyword search query to input into NCBI In order to find all the available proteins accession numbers for a specific family?

For example, how would I find all the reviewed accession numbers for all currently available human ALDP (AdrenoLeukoDystrophy Protein) family of proteins that are in NCBI?

Also, is there a specific database within NCBI that is optimized for my task?

Thanks.

# Answer
> 2 votes
A simple approach is to search for the gene name in NCBI Gene portal but this can miss a few cases and often include genes that are unrelated because the gene symbol is used by other unrelated genes as well.

For a cleaner dataset, you can use the NCBI Orthologs as follows. Search for `human ALDP` in the NCBI homepage with `All Databases` selected in the dropdown menu as shown below:

After you click the 'Search' button, you will be directed to a result page that looks like this:

Click on the 'Orthologs' button and you will be directed to a page with the list of proteins that are similar to the one you have started from. Note, these may not be orthologs in the truest sense; see this description about how they are calculated.

You can search for the taxonomic group that you are interested in, pick and choose the organisms of your interest and download the table in CSV format. It will have the NCBI GeneID, RNA and protein accessions.

# Answer
> 0 votes
*This is more of a comment, but due to the lengthiness I am posting it as an answer, which it isn't.*

If by protein family you mean more than the 3-5 similar protein in bony fish and tetrapods resulting from the x4 genome duplication after the split with lampreys, the solution by @vkkodali is perfect for you —and I think that is the case. Doing a Blast of the gene with `human` and `Petromyzon marinus` as an outgroup also works well especially if you want to infer some relationship between them. The huge ABC transporter family is one of the exceptions, so this trick won't help as it duplicate many many times across the tree of life.

If you meant *domain fold* family, *i.e.* in the structural sense, like TIM barrels, Rossmann folds, β-propellers (WD40 repeat), SH3 and so forth, you are after PFam families, for which Uniprot is a much better database than NCBI, albeit trickier to parse programmatically: uniprot: P33897 gives Pfam: PF00005.

---
Tags: phylogenetics, literature-search
---


# Question
Title: What kind of analysis is practically done on GSE data files?
I have a GSE data file in csv file format containing fields such as: ID, adj.P.Val, P.Value, t, B, logFC, Gene.symbol, Gene.title. In which adj.P.Val, P.Value, t, B, logFC fields being numeric. What are the factors I need to consider if I want to cluster the data only on logFC using K-Means clustering algorithm ? And first of all is it feasible to perform clustering on GSE data files ? If yes, what should be the approach ? If not, what different kinds of analysis can be performed on such kind of datasets ?

Following is snapshot of the data:

# Answer
Have you ever used `GEO2R`? You can use `GEO2R` to analyze the data dividing it in subgroups (clusters) that interest you.
> 1 votes

# Answer
You have multiple questions:

First of all, here is the description of `GSE files` or `GEO Series`:

> **Series records are supplied by submitters**
>
> A Series record links together a group of related Samples and provides a focal point and description of the whole study. Series records may also contain tables describing extracted data, summary conclusions, or analyses. Each Series record is assigned a unique and stable GEO accession number (GSExxx).

Apparently the authors have provided their (differential expression) analysis results (p-values, adjusted p-values, logFC, ...). And the table you are sharing in your question does not look like a typical GSE profile.

On your question on the feasibility of clustering based on GSE data: that would depend on the features (columns) you are provided with, you should have an understanding of each feature that you have. In your case, you can cluster your genes based on p-value, adjusted p-value and logFC, but would that make sense?

Regarding your question using k-means on only one variable, logFC: This link points to a related discussion. Short answer is "technically yes".

Your last question depends on the "research question" you have, if you are asking about the table in your question: pathway analysis, gene set enrichment analysis, checking for co-expressed modules, ...
> 1 votes

---
Tags: bioconductor, python, clustering, machine-learning
---


# Question
Title: Which correlation method to compute the correlation score between different clusters of Sc-RNAseq data?
I'm analyzing single cell rna-seq data and trying to compute the correlation score between different clusters. Wondering how to choose the correlation method("pearson" (default), "kendall", or "spearman") of the Cor function?

I'm a new comer in bioinformatics and know little about these arithmetics. Can anyone give me some advice or show me the code to decide which method to choose for my data matrix?

Below is my gene expression matrix: rows are gene names,and columns are clusters.I use the default Peasrson method to do analysis,but I really don't know if it is right.

```
CorOb.cor.exp <- as.data.frame(cor(CorOb.av.exp))
CorOb.cor.exp$x <- rownames(CorOb.cor.exp)
neupeptide.cor.df <- tidyr::gather(data = CorOb.cor.exp, y, correlation,CorOb.cor.exp$x)
a <- ggplot(neupeptide.cor.df, aes(x, y, fill = correlation)) + geom_raster()+scale_fill_gradientn(colors = c("blue", "white", "red"))+
  theme(axis.text.x = element_text(angle = 45, vjust = 1,
                                   size = 12, hjust = 1))+
  coord_fixed()

```

a

# Answer
Disclaimer: This is not a complete answer but more a comment, posting as an answer because of character restrictions for comments. Also, I am not a statistican but a molecular biologist who did quite some reading/courses on stats.

Two scRNA-seq tools, namely `SingleR` and `RCA` (Reference Component Analysis, Li et. al, 2017) compute correlations of bulk transcriptomic data to single cells; for this they use `Spearman correlation` and `Pearson correlation` respectively (The two tools use the same bulk transcriptome reference). I am using these two tools for different purposes so I did not really compare their results, however, single cell clusters that I compute using the `reference components` (calcuated with `RCA`) is in good agreement with single cells labels computed with `RCA` (adjusted rand index of 0.81).

In addition to that, please take a look at scMatch: a single-cell gene expression profile annotation tool using reference datasets, they claim (scMatch computes bulk - single cell correlations like SingleR and RCA):

> From this comparison we found Spearman’s correlation coefficient outperformed Pearson’s correlation for all six libraries.
> 1 votes

# Answer
This question does not have a definite answer. Rather, the answer depends heavily on the assumptions you wish to make regarding your data.

Given a pair of random variables $x$ and $y$, Pearson's rho has a value of 1 when there is a linear relationship between $x$; i.e. $y = x\beta + c$, where $c$ is a constant, and there is no noise. Note that the slope, $\beta$, and intercept, $c$, of this relationship do not affect the correlation coefficient. With Spearman's rho, we first compute the rank of each element in $x$ and $y$. We then calculate the Pearson correlation between these ranks.

Kendall's tau correlation coefficient also deals with the ranks of $x$ and $y$, but is slightly more complex. The most relevant difference in this case is likely that tau is better in the case of ties, and may be more robust for low sample sizes. In most cases, Kendall's and Spearman's are likely to give very similar results

If you wish to assess only linear relationships between expression profiles, Pearson's is the most appropriate. If you wish to allow for non-linear effects, Spearman or Kendall is most appropriate (Kendall if the data contains ties).
> 1 votes

# Answer
Pearson correlation evaluates two continuous variables. Spearman correlation is often used to evaluate ordinal variables. Kendall correlation is a nonparametric test for two continuous variables.
> 0 votes

---
Tags: r, rna-seq, scrnaseq, seurat, single-cell
---


# Question
Title: multiBamSummary output chromosome region annotate issue
I was looking how to take readcount from ATAC seq data , this post which suggested to use bedtools and featurecounts to take out readcounts .But multiBamSummary seems to be pretty straight forward. Now I have my count file but only issue is there are no gene names. So now is it possible to annotate these chromosome regions? if so how ?

Any suggestion or help would be highly appreciated

# Answer
> 1 votes
Assuming you used a BED file as input then the output is in the same order, so just merge the two files together.

---
Tags: deeptools
---


# Question
Title: Interpreting a kernel density function plot of gene expression
I need help interpreting a kernal density plot and a histogram of this data. I have provided some context below let me know if you need more information

A co-worker that I am collaborating with has a list of genes ranked by expression fold change. He says most of the genes in one dataset are noise with respect to his analysis so he is trying to filter this out. He has sent me a document of his methods for generating a new list of genes which he think are relevant. he has run a permutation test on the ranked genes I have provided a quote here:

"To determine which genes contributed more strongly in real data, as compared by permuted data, the rankings achieved for each gene in real data were compared to the shuffled rankings. For each gene, and each permutation type, a p-value was calculated by counting how many times that gene was ranked higher in shuffled data than in real data, and dividing it by the number of comparisons (16 * 256 = 4096). To achieve p\<0.05 a gene’s real ranking would have to be higher than 95% of its shuffled rankings."

Following this he provides kernel density plots of the genes under different rankings.

"Below are kernel density estimation (kde) plots, overlaid on histograms, of the rankings of two genes randomly selected from the few genes with p\<0.05 vs all three permutations."

"Below are kde plots of the rankings of two genes randomly selected from the many genes with p\>0.05 vs all three permutations."

What am I looking at? What does this tell me about the genes ranking statistics.

# Answer
There are two topics here to be discussed. Gene Set Enrichment Analysis and a kernel Density plot.

## GSEA

The Wikipedia page about Gene Set Enrichment Analysis covers well the principle and the tools that can be used to calculate it. To cite it:

> The general steps are summarized below:
>
> Calculate the enrichment score (ES) that represents the amount to which the genes in the set are over-represented at either the top or bottom of the list. This score is a Kolmogorov–Smirnov-like statistic \[a ranking based statistic\].
>
> Estimate the statistical significance of the ES. This calculation is done by a phenotypic-based permutation test in order to produce a null distribution for the ES. The P value is determined by comparison to the null distribution.
>
> Calculating significance this way tests for the dependence of the gene set on the diagnostic/phenotypic labels

Basically, you want to test if the ranking is random or not. Therefore you make a random dataset and test whether your set is like that. One assumption underlying this is that the output of the gene whatever it may be is *not* normally distributed (parametric, follows a bell curve), hence the use of ranking, which removes this (you lose signal but you lose biases). I am assuming the non-parametric nature of your ouput values was tested or is a known non-parametric value, such as RPKM or log2 fold change (by virtue of being a ratio). In other words, the magnitude difference between your samples is lost, but their order is not.

There are several tools, some online, that do this enrichment testing for you depending on the task, some even work with gene ontologies.

## KDE

A kernel density plot is a like a histogram, but smoothed, albeit not in a moving average way. In a histogram you divide your x-axis into bins: it is discreet and you get an integer count per bin. In a kernel density plot the data is fitted to "probability density function", an equation which given x will give y, where the integral of the curve is one, hence why y is density not counts. The x axis will be your bins-like value, in this case rank. The binning window metric for a KS function is called bandwidth.

See Matlab documentation for a fairly non-opaque explanation.

However, if you want to compare histograms, I strongly suggest violin plots, not because they are strongly in vogue, but because you can compare histograms without worrying about overlaying them.
> 2 votes

---
Tags: statistics, gene-expression, gene-ontology
---


title: CRISPR screening and systems biology

text:
# Question
Title: CRISPR screening and systems biology
I am not sure how to tackle this one, so I’ll explain the general idea I have in mind.

Given results of a CRISPR knockout screening experiment (like in this paper) with two types of assays: control vs. treated with 5 time points per each of the assays.

Is there a way to use the screening data to perform systems biology network analysis? Something which is related to signal transduction, transcriptional regulatory, or metabolic networks?

Maybe something which uses Bayesian inference? Or something similar to differentially expressed representation of the CRISPR screening? Alternatively, a GWAS analysis similar to this paper?

# Answer
What do you call screening results? Do you have a kind of array with protein or RNA data? If so, you can use several open-source tools available on the internet, like Enrichr or network analyst to better understand the pathways affected by your treatment.
> 1 votes

# Answer
To investigate genes through a network based approach you can use several tools. I don't know which level of programming you have or which tool you prefer but here is a brief overview of the possibilities you have.

You can look into established databases containing information about interactions like, for instance, STRING. Or, you can build networks using the wgcna method with R.

Then, once you have your network you can analyse its properties and here you also have several tools available such as standalone, java-base software like cytoscape, R libraries such as igraph, python libraries such as networkX online tools like centiserver...

Also you can perform gene set enrichment analsysis, which means adding a further layer of biological information to you set of genes like metabolic pathways (from KEGG, for example) and molecular functions (from GeneOntology, for instance) (gsea) as here, here, and here.

Here there is also a nice tutorial about network analysis with R.

---
> 1 votes

---
Tags: genomics, gwas, crispr, systems-biology
---


# Question
Title: Converting Gene Symbol to Ensembl ID in R
I'm trying to convert ~20,000 different human gene symbols to ensembl IDs. I've been trying to use biomaRt to do this, but continue getting the following error

```
getBM( attributes=c("ensembl_gene_id") , filters= "mgi_symbol"  ,mart=ensembl)
Error in martCheck(mart) :
  No dataset selected, please select a dataset first.  You can see the available datasets by using the listDatasets function see ?listDatasets for more information.  Then you should create the Mart object by using the useMart function.  See ?useMart for more information

```

Is there a different format I should be using for this? I tried using BioConductor's description of how to use getBM but am a R novice and confused as to how to set this up to convert gene symbols into ensemble IDs.

# Answer
You can use the standard annotation package for humans (you can also use biomaRt, but it can be more confusing, see below):

```
library("org.Hs.eg.db") # remember to install it if you don't have it already
symbols <- mapIds(org.Hs.eg.db, keys = ensemblsIDS, keytype = "ENSEMBL", column="SYMBOL")

```

The `ensemblsIDS` should be a character vector: `c("ENSG001", "ENSG002", ..)`.

However, if your primary key is an ENSEMBL id it is preferred to use a ENSEMBL package centred around the ENSEMBL id :

```
library("EnsDb.Hsapiens.v86")
mapIds <- mapIds(EnsDb.Hsapiens.v86, keys = ensemblsIDS, keytype = "GENEID", columns = c("SYMBOL")

```

The solution with biomaRt:

```
library("biomaRt")
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
getBM(attributes='hgnc_symbol',
      filters = 'ensembl_gene_id',
      values = ensemblsIDS,
      mart = ensembl)

```
> 12 votes

# Answer
I tried several R packages (`mygene`, `org.Hs.eg.db`, `biomaRt`, `EnsDb.Hsapiens.v79`) to convert Ensembl.gene to gene.symbol, and found that the `EnsDb.Hsapiens.v79` package / gene database provides the best conversion quality (in terms of being able to convert most of Ensembl.gene to gene.symbol).

```
# Install the package if you have not installed by running this command:
BiocManager::install("EnsDb.Hsapiens.v79")

library(EnsDb.Hsapiens.v79)

# 1. Convert from ensembl.gene to gene.symbol
ensembl.genes <- c("ENSG00000150676", "ENSG00000099308", "ENSG00000142676", "ENSG00000180776", "ENSG00000108848", "ENSG00000277370", "ENSG00000103811", "ENSG00000101473")

geneIDs1 <- ensembldb::select(EnsDb.Hsapiens.v79, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))

# 2. Convert from gene.symbol to ensembl.gene
geneSymbols <-  c('DDX26B','CCDC83',  'MAST3', 'RPL11', 'ZDHHC20',  'LUC7L3',  'SNORD49A',  'CTSH', 'ACOT8')

geneIDs2 <- ensembldb::select(EnsDb.Hsapiens.v79, keys= geneSymbols, keytype = "SYMBOL", columns = c("SYMBOL","GENEID"))

```

Other available R packages / gene database for the conversion can refer to this GitHub page.
> 6 votes

---
Tags: r, gene, ensembl, biomart
---


title: How do you map nanopore fast5 files?

text:
# Question
Title: How do you map nanopore fast5 files?
Is there a best practices for mapping Oxford Nanopore files to a reference?

Is there a tool that can take a tarball of fast5 files and map them directly or do they need to be converted to fastq first?

# Answer
> 2 votes
Fast5 files need to be converted to fastq first prior to mapping. There are currently no tools available to map the entire raw signal directly to a reference, although a few tools are available for identifying short sequences within the raw signal.

---
Tags: read-mapping, nanopore, fast5
---


title: How to do the reverse of bcf merge?

text:
# Question
Title: How to do the reverse of bcf merge?
How can one do the reverse of bcftools merge? To get one sample back out of the merged one. Without holding reference sites.

# Answer
> 1 votes
`bcftools` has the `view` command which can do this. Create a text file (`samples`) which has the names of the samples you want, one per line, and then:

```
bcftools view -S samples merged.vcf.gz > new.vcf

```

Or, if you're running a shell that supports it, you can use input redirection for a single sample name:

```
bcftools view -S <(echo "sample name") merged.vcf.gz > new.vcf

```

---
Tags: bcftools
---


title: SRA Toolkit and lebanese data

text:
# Question
Title: SRA Toolkit and lebanese data
I am trying to extract data from this: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR6245218

I installed SRA Toolkit, downloaded the SRR6245218 file and executed this:

```
fastq-dump -X 5 SRR6245218 --fasta -O ./

```

The result is a fasta file 900 bytes long, which does not contain the mitochondrial data that I need. What am I doing wrong?

# Answer
Your command is fine. It is working as expected. Did you intend to include `-X 5` in your command? That restricts the number of spots to 5. If you drop it, you will download the entire set of 2563898 spots.
> 1 votes

---
Tags: fasta, phylogenetics, sratoolkit
---


title: How to use -csvStats option?

text:
# Question
Title: How to use -csvStats option?
I'm trying to get CSV files from snpEff -csvStats option and got stuck. I tried the code below and only got the vcf file with no contents. How can I use -csvStats option properly?

```
snpEff -Xmx4g -csvStats -v hg19 08_FILTER/NA12878.Filtered.Variants.vcf > test.ann.vcf

    00:00:00  SnpEff version SnpEff 4.3t (build 2017-11-24 10:18), by Pablo Cingolani
    00:00:00  Command: 'ann'
    00:00:00  Reading configuration file 'snpEff.config'. Genome: '08_FILTER/NA12878.Filtered.Variants.vcf'
    00:00:00  Reading config file: /home/ubuntu/01_NA12878/snpEff.config
    00:00:00  Reading config file: /home/ubuntu/miniconda2/share/snpeff-4.3.1t-2/snpEff.config
    java.lang.RuntimeException: Property: '08_FILTER/NA12878.Filtered.Variants.vcf.genome' not found
      at org.snpeff.interval.Genome.<init>(Genome.java:106)
      at org.snpeff.snpEffect.Config.readGenomeConfig(Config.java:681)
      at org.snpeff.snpEffect.Config.readConfig(Config.java:649)
      at org.snpeff.snpEffect.Config.init(Config.java:480)
      at org.snpeff.snpEffect.Config.<init>(Config.java:117)
      at org.snpeff.SnpEff.loadConfig(SnpEff.java:451)
      at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:1000)
      at org.snpeff.snpEffect.commandLine.SnpEffCmdEff.run(SnpEffCmdEff.java:984)
      at org.snpeff.SnpEff.run(SnpEff.java:1183)
      at org.snpeff.SnpEff.main(SnpEff.java:162)
    00:00:00  Logging

00:00:01  Done.

 snpEff -Xmx4g -v hg19 -csvStats 08_FILTER/NA12878.Filtered.Variants.vcf > tes2.ann.vcf

00:00:25  Genome stats :
#-----------------------------------------------
# Genome name                : 'Homo_sapiens (USCS)'
# Genome version             : 'hg19'
# Genome ID                  : 'hg19[0]'
#     'M' 16571   Vertebrate_Mitochondrial
#     'Un_gl000226'   15008   Standard
#     '18_gl000207_random'    4262    Standard
#     'MT'    1   Vertebrate_Mitochondrial
#-----------------------------------------------

00:00:27  Predicting variants

```

Here is the snpEff.config file. I installed the snpEff at /miniconda2/share .

```
# Databases are stored here
# E.g.: Information for 'hg19' is stored in data.dir/hg19/
#
# You can use tilde ('~') as first character to refer to your home directory.
# Also, a non-absolute path will be relative to config's file dir
#
#---
data.dir = ./data/

#---
# Database repository: A URL to the server where you can download databases (command: 'snpEff download dbName')
#---
database.repository = http://downloads.sourceforge.net/project/snpeff/databases

#---
# Latest version numbers. Check here if there is an update.
#---
versions.url = http://snpeff.sourceforge.net/versions.txt

#-------------------------------------------------------------------------------
# Third party databases
#-------------------------------------------------------------------------------

# Databases for human GRCh38 (hg38)

database.clinvar.GRCh38                 = ./db/GRCh38/clinvar/clinvar-latest.vcf.gz
database.repository.clinvar             = ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar-latest.vcf.gz

database.dbsnp.GRCh38                   = ./db/GRCh38/dbSnp/dbSnp.vcf.gz
database.repository.dbsnp.GRCH38        = ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/00-All.vcf.gz

database.dbnsfp.GRCH38                  = ./db/GRCh38/dbNSFP/dbNSFP.txt.gz

# Databases for human GRCh37 (hg19)
database.clinvar.GRCh37                 = ./db/GRCh37/clinvar/clinvar-latest.vcf.gz
database.repository.clinvar.GRCh37      = ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh37/clinvar-latest.vcf.gz

database.dbsnp.GRCh37                   = ./db/GRCh37/dbSnp/
database.repository.dbsnp.GRCh37        = ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b141_GRCh37p13/VCF/

database.gwascatalog.GRCh37             = ./db/GRCh37/gwasCatalog/gwascatalog.txt
database.repository.gwascatalog.GRCh37  = http://www.genome.gov/admin/gwascatalog.txt

database.dbnsfp.GRCh37                  = ./db/GRCh37/dbNSFP/dbNSFP.txt.gz

#-------------------------------------------------------------------------------
# Loss of function (LOF)
#-------------------------------------------------------------------------------

# It is assumed that even with a protein coding change at the
# last 5% of the protein, the protein could still be functional.
lof.ignoreProteinCodingAfter  : 0.95

#  It is assumed that even with a protein coding change at the
#  first 5% of the protein:
#       "..suggesting some disrupted transcripts are
#       rescued by transcriptional reinitiation at an
#       alternative start codon."
lof.ignoreProteinCodingBefore : 0.05

# Larger deletions removing either the first exon or more than
# 50% of the protein-coding sequence of the affected transcript
lof.deleteProteinCodingBases : 0.50

#-------------------------------------------------------------------------------
# Codon tables
#
# Format:   It's a comma separated "codon/aminoAcid[+*]" list
#           Where 'codon' is in uppper case, aminoAcid is a one letter
#           code, '+' denotes start codon and '*' denotes stop codon.
#
# References:   http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi
#               ftp://ftp.ncbi.nih.gov/entrez/misc/data/gc.prt
#-------------------------------------------------------------------------------

```

# Answer
Had the same issue, but found that I need to provide a filename to use the csvstats option

```
snpEff -Xmx4g -csvStats stats.csv -v hg19 08_FILTER/NA12878.Filtered.Variants.vcf > test.ann.vcf

```
> 0 votes

---
Tags: variation, gatk, snpeff
---


# Question
Title: How to quickly and robustly convert between mmCIF and PDB?
There is already a question on PDB/CIF to MMTF, however what is a robust way to programmatically go between PDB and CIF files?

For example I can use a python script from this gist that relies on Biopython.

```
id="4ckh"
wget http://www.ebi.ac.uk/pdbe/static/entry/download/$id-assembly-1.cif.gz
gunzip $id-assembly-1.cif.gz
python cif2pdb.py $id-assembly-1.cif $id.pdb

```

This process yields a pdb structure as shown below.

The issues are:

* Requires python/Biopython to be installed.
* The output isn't perfectly clean.
* The gist is relatively undeveloped and hasn't got much evidence of testing.

# Answer
> 5 votes
Converting a PDB file to some version of an mmCIF file is always possible provided the PDB file doesn't have any "problems". However grouping identical molecules into the same entity is non-trivial so the Biopython and BioJulia mmCIF writers (which I wrote) do not try and do this, simply assigning new entites for each molecule. These files will also be missing other things like the metadata parts of the mmCIF dictionary not involving the atoms.

In Biopython:

```
from Bio.PDB import PDBParser, MMCIFIO
p = PDBParser()
struc = p.get_structure("", "file.pdb")
io = MMCIFIO()
io.set_structure(struc)
io.save("file.cif")

```

In BioJulia:

```
using BioStructures
struc = read("file.pdb", PDB)
writemmcif("file.cif", struc)

```

Going the other way is, however, not always possible. You can run out of atom serials or chain IDs to use in the output PDB file (this is one reason for the switch to mmCIF files in general). The gist you have linked is a way to get round this.

If this is not a problem for your structure you can invert the above snippets to get a decent PDB file - see the relevant docs. Note that in both cases it will use the `auth_seq_id` and `auth_asym_id` fields to get the residue number and chain ID rather than the relabelled `label` versions. The test suites of Biopython and BioJulia do test these conversions pretty thoroughly.

There is also pdb-tools if you do not want to install Biopython, though both involve a simple `pip install` so it's not much difference.

For some speed benchmarks of various software see here and here.

---
Tags: file-formats, protein-structure
---


# Question
Title: How can I run a command for multiple files?
I have a bunch of `.vcf` files in a folder and I want to run command below on all of them but doing that one by one manually is really painful. I am seeking for a way in terminal to do that simultaneously. The output would be a `.txt` file.

```
grep -v "##" /file.vcf | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > /file.txt

```

Can you help please?

# Answer
Assuming all your vcfs are in the same folder:

```
for file in *.vcf; do
    grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt
done

```
> 5 votes

# Answer
If additional CPUs are available:

```
CPU_COUNT=5
for file in *.vcf; do
     (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1
     [[ $((count%CPU_COUNT)) -eq 0 ]] && wait
done

```

one-liner:

```
CPU_COUNT=5; for file in *.vcf; do (grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file%.vcf}.txt) & let count+=1; [[ $((count%CPU_COUNT)) -eq 0 ]] && wait; done

```
> 2 votes

---
Tags: bash, linux, loop
---


# Question
Title: how to solve “dim(X) must have a positive length” at running ComBat function in R
I used `ComBat()` for batch effect correction in my expression data. basically, that function inputs are expression data, Batch covariate, and Model matrix for the outcome of interest and other covariates besides batch. So, I prepared all inputs based on the instruction of `ComBat()` function in 'sva' package. "Pheno\_LMS" is my clinical data inclus batchId of my samples and "TCGA\_expr\_log" is my expression data. my code is below:

Pheno\_LMS structure:

```
Sample        batchId   age  . . .
GSM71019.CEL    396     63
GSM71020.CEL    396     58
GSM71021.CEL    410     85
GSM71020.CEL    411     58
GSM71021.CEL    410     40
.
.
.
dim(Pheno_LMS)
[1] 74 37

```

TCGA\_expr\_log structure:

```
                   GSM71019.CEL 1      GSM71019.CEL 2        GSM71019.CEL 3  . . .
    Gene id
  ENSG00000000003     10.115170           8.628044                  8.779235
  ENSG00000000005     5.345168            5.063598                   5.113116
  ENSG00000000009     6.348024            6.663625                   6.465892
    .
    .
    .
dim(TCGA_expr_log)
[1] 42817    74

batch = Pheno_LMS$batchId
modcombat = model.matrix(~1, data=Pheno_LMS)
TCGA_expr_Co <-  ComBat(as.matrix(TCGA_expr_log),batch = batch,mod = modcombat, par.prior = TRUE,mean.only = TRUE)

```

but when I run that function, I got below error:

```
Error in apply(dat[, batch == batch_level], 1, function(x) { :
  dim(X) must have a positive length

```

Now, I need anybody comments about my problem. thanks

P.S: when i call batch I see:

```
> batch
 [1] 13 13 13 14 14 14 14 14 14 14 14 6  6  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
[29] 7  8  7  8  11 11 11 14 12 1  6  4  2  3  3  3  3  6  6  1  3  9  8  6  11 7  3  5
[57] 4  5  5  4  4  11 8  8  12 12 12 12 12 10 10 10 11 13
Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14

```

# Answer
`batch = Pheno_LMS$batchId` should be `batch = Pheno_LMS$BatchId`, if you fix that then at least the code you showed works just fine.
> 1 votes

# Answer
as suggested in the comment you should verify if the command:

```
dat[, batch == batch_level]

```

returns an empty object or not. Apparently, *batch* and *batch\_level* do not have anyting that match and the function tries to run over an empty dataframe.
> 0 votes

---
Tags: r, batch-effects, combat, sva
---


title: Add new tool to galaxy

text:
# Question
Title: Add new tool to galaxy
I am trying to include a new tool to the galaxy main menu, following these instructions.
`tool_conf.xml` doesn't exists only `tool_conf.xml.main` and `tool_conf.xml.sample` and none of them seem to be working. I make changes and nothing happens when I run galaxy.

What am I missing?

Thank you.

# Answer
> 1 votes
The .main and .sample file files usually come with default installations and may have been removed by whoever setup this Galaxy instance. You can look at the originals online and just base your own file on those: https://github.com/galaxyproject/galaxy/blob/dev/lib/galaxy/config/sample/tool\_conf.xml.sample

---
Tags: ngs, galaxy
---


# Question
Title: Can I use a regular liner regression model when I'm working with DNA methylation data?
I'm new on working with genetics data and I'm just wondering if I can use the `lm` function when I'm building my models or do I have to use `lmFit` function in limma package to build the models and what is the difference between the two packages in terms of performing linear regression models

# Answer
You can use either, but `lmFit` has the benefit of returning an object that can be used with `eBayes()` so you can pool information across genes/probes/whatever. `lm()` is a base R function applicable basically everywhere. `lmFit()` if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with it since sample sizes tend to be small so everyone want's to pool information for more statistical power.

BTW, if you have methylation data then you'll want to fit data in logit space.
> 2 votes

---
Tags: limma, linear-regression
---


# Question
Title: Remove repeated ALT allele from the REF field of a vcf-like file
I have 200 `.txt` files on a Linux system. In each file, I have the reference allele in the `REF` column and INDELs in the `ALT` column. Let's say I have `A` in the REF column and `AC` in the ALT. I want to remove the `A` from the ALT column because that is a repetition of the `A` reference allele. For example, if I have `C` in the REF column and `CCAT` in the ALT column, I only want CAT to remain in the ALT column.

This is what my files look like

```
CHROM   REF ALT
chr1    C   CT
chr1    C   CT
chr1    C   CA
chr1    C   CT

```

Can someone help me please?

Solution from @Terdon

This is my original file

```
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   GATAAC
chr4    30280292    .   C   CACGTAT
chr6    166253422   .   C   CCTTTCT
chr11   21931340    .   A   AACACACGC

```

This is what happened after running your kindly suggestion

```
[fi1d18@cyan01 indel]$ awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,4); }}1' indels.vcf > out.vcf

#CHROM  POS ID  REF ALT
chr4    5708160 .   G   AAC
chr4    30280292    .   C   GTAT
chr6    166253422   .   C   TTCT

```

I think more than one allele being removed from ALT column, is not it?

# Answer
I am making a couple of assumptions here:

1. the delimiter is tab
2. the reference allele is always present in the ALT column
3. the reference allele in the ALT column is always in the first position

With those in mind, you can do this using `awk` as follows:

```
awk 'BEGIN{FS="\t";OFS="\t"}{$3=substr($3,2)}{print $0}'

```

Update 1: Using the new VCF style file that @Angel provided, here's how I would do it. This is not significantly different from what @terdon proposed below. The code below will work even if the reference allele is \>1 nt.

```
$ cat test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       GATAAC
chr4    30280292        .       C       CACGTAT
chr6    166253422       .       C       CCTTTCT
chr11   21931340        .       A       AACACACGC

$ awk 'BEGIN{FS="\t";OFS="\t"}{l=length($4); if ($4==substr($5,1,l)) {$5=substr($5,l+1); print $0} else {print $0}}' test.vcf
#CHROM  POS     ID      REF     ALT
chr4    5708160 .       G       ATAAC
chr4    30280292        .       C       ACGTAT
chr6    166253422       .       C       CTTTCT
chr11   21931340        .       A       ACACACGC

```
> 4 votes

# Answer
Another awk approach. This should work on a normal VCF file (it will print the headers) as well as your example and only changes the ALT field if its first character is the same as the 2nd field:

```
$ awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,2); }}1' ff.vcf
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   ATAAC
chr4    30280292    .   C   ACGTAT
chr6    166253422   .   C   CTTTCT
chr11   21931340    .   A   ACACACGC

```

To run this on all of your files, assuming they are the only files in the current directory with a `.txt` extension, you can do:

```
for file in *txt; do
    awk -F"\t" -vOFS="\t" '{ if(substr($5,1,1) == $4){$5=substr($5,2); }}1' "$file" > "fixed.$file"
done

```

This will loop over all `*.txt` files in the current directory. For each `foo.txt`, it will create `fixed.foo.txt` with the modified data.
> 2 votes

---
Tags: vcf, text-processing
---


title: WGCNA module preservation doubt

text:
# Question
Title: WGCNA module preservation doubt
The WGCNA analysis page48 of that link.

In this analysis they are comparing the Preservation of modules between human and chimpanzee brain networks where i see the number of human modules and chimp modules are exact same in number as well as the color they are labeled.

So is it necessary that two different condition where, Im comparing need to have equal number of modules to check the module preservation? Because the power would be different,then the downstream output would be different for two set of condition ,let say if Im trying to see what are the module conserved in my Control vs a disease condition ?

I have 4 different condition such as CONTROL,BPD,SCZ and MDD all\[micro-array expression data\] brain disorder with samples number of 50,54,54,48.

So what im trying to see is what are the modules are getting preseved in each diseased state when they are compared to Control condition.

What i have done is I ran WGCNA on each condition separately and found what are the modules are there in different condition. So here I see the number of modules that are observed in each individual condition varies.

Now if I do module preservation is it necessary to have similar number of modules from each condition ? Or its fine to have different module number from different condition and can still run the module preservation.

These are gene expression network. Modules are obtained as explained in the WGCNA tutorial.

**Do you mean "conservation" or "preservation"?** I mean module preservation

Suggestion would be highly appreciated

# Answer
The networks you want to compare do not have to have the same number of modules. (The human and chimp networks just happen to have the same number of modules.) In fact, you only need modules in the reference network; modules from the test network are not used at all and need not be defined. You may want to read the original article describing module preservation, it would help you understand how the method works, what it needs and what it does not.
> 1 votes

---
Tags: wgcna
---


# Question
Title: Illustrate a 3D visualisation of the three main PCs using plot3d() package in R?
I have a dataset (view/download) here. I need to run PCA on this data set and need to illustrate a 3D visualisation of the three main PCs using plot3d() package.

I have looked for more information before coming here on which package to use, some people use princomp and some us prcomp (in university we have personally been using prcomp). The confusing part is were to move forward with this?

The full question is "Illustrate a 3D visualisation of the three main PCs using plot3d() package in R. Use the following colours for samples (data points) belonging to the various immune subgroups/subtypes (C1: red, C2: yellow, C3: green, C4: cyan, C5: blue, C6: purple)"

After reading in the CSV and loading rgl library

In the current state I have removed the rows with the subgroups: see code below

```
# Remove last row from dataset
dataMinusLastRow = data[1:(nrow(data)-1),]

# Convert dataset to numeric and apply
coltemp <- colnames(dataMinusLastRow)
pcaData = t(apply(dataMinusLastRow, 1, as.numeric))
colnames(pcaData) <- coltemp

# PCA
pc = prcomp(t(pcaData), scale=T, center=T)
summary(pc) # From this we can see that PC184 accounts for 95% of cumulative proportion

# Getting the cumulative proportion in numeric dataset
vars <- apply(pc$x, 2, var)
props <- vars / sum(vars)
pc.cs <- cumsum(props) # Cumulative Proportion
pc.index<-min(which(pc.cs>0.95)) # Index of which cumulative proportion > 95%

```

I think the main part which confuses me is that it asks for the three main PC's which I am not sure what they are and have had not much support from university with information regarding this.

Thanks to anyone who can help in any shape/form!

# Answer
Regarding your question on how to decide which PCs to choose for plotting, `prcomp` is ordering PCs by their proportions in the variance. For example, using the dataset `iris`:

```
iris_pc <- prcomp(iris[,1:4], scale  = TRUE, center = TRUE)

```

and you get

```
> summary(iris_pc)
Importance of components:
                          PC1    PC2     PC3     PC4
Standard deviation     1.7084 0.9560 0.38309 0.14393
Proportion of Variance 0.7296 0.2285 0.03669 0.00518
Cumulative Proportion  0.7296 0.9581 0.99482 1.00000

```

BTW, using the `princomp` instead of `prcomp` gives similar results:

```
iris_pc2 <- princomp(iris[,1:4], cor=TRUE, scores=TRUE)

```

```
> summary(iris_pc)
Importance of components:
                          Comp.1    Comp.2     Comp.3      Comp.4
Standard deviation     1.7083611 0.9560494 0.38308860 0.143926497
Proportion of Variance 0.7296245 0.2285076 0.03668922 0.005178709
Cumulative Proportion  0.7296245 0.9581321 0.99482129 1.000000000

```

So, this means if you want to represent the main 3 PCs, you have to pick the first three PCs.

However, based on your dataset, if you take a look at the output of `prcomp`, you will see:

```
pc = prcomp(t(pcaData), scale=T, center=T)

```

```
> summary(pc)
Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6     PC7    PC8     PC9    PC10   PC11    PC12
Standard deviation     9.3653 7.9806 5.46230 5.02379 3.95681 3.36495 3.21719 3.1042 2.98449 2.73437 2.6529 2.56841
Proportion of Variance 0.1993 0.1447 0.06781 0.05736 0.03558 0.02573 0.02352 0.0219 0.02024 0.01699 0.0160 0.01499

```

You can see that the first three PCs count for only 41% of the total variance. This can be explain by the diversity of your samples (coming from different cancers for example). So it means that your graphic representation won't be perfect to separate clusters.

Here, the graphic representation of the cumulative proportion of variances accross PC:

```
plot(pc, type = "l")

```

Next for the plotting and color attribution, I get a lot of trouble to extract the values from your `Subgroups` row.... for some reasons, I was not able to convert them into a character vector... weird ... anyway, I re-defined groups based on the code you put in patients ids:

```
patients = as.character(colnames(data))
grp = c(rep("C1",length(grep(".C1_",patients))),
        rep("C2",length(grep(".C2_",patients))),
        rep("C3",length(grep(".C3_",patients))),
        rep("C4",length(grep(".C4_",patients))),
        rep("C5",length(grep(".C5_",patients))),
        rep("C6",length(grep(".C6_",patients))))

```

For plotting results of PCA, you can use plot3d from `rgl` package for example. It gives you an interactive plot.

```
library(rgl)
colors = c("red","yellow", "green","cyan", "blue","purple")
plot3d(pc$x[,1:3], col =colors[as.factor(grp)], main = "pc / prcomp")
legend3d("topright",legend = unique(grp),col = colors, pch = 16, cex = 2)

```

And here the snapshot of your plot:

Hope that it helps you to figure it out the solution of your problem
> 1 votes

---
Tags: r, pca
---


title: Modifying .vcf files

text:
# Question
Title: Modifying .vcf files
I have 200 .vcf files like this on a Linux system:

```
#CHROM  POS ID  REF ALT
chr4    5708160 .   G   GATAAC
chr4    30280292    .   C   CACGTAT
chr6    166253422   .   C   CCTTTCT
chr11   21931340    .   A   AACACACGC
chr12   118404422   .   G   GCACACA
chr15   22672047    .   T   TAAGTAGATTTAAG
chr16   32656293    .   ATGATT  A

```

I want to modify my .vcf files in this way

For instance if the reference allele in `REF column` is `TCTC` and alteration allele in `ALT column` is `T`, because T is the first common nucleotide in both reference and alteration so we remove `T` from `ALT column` so REF would remain as `TCTC` and ALT would be `-` . In total I mean removing whatever in ALT which is common with the first allele in REF.

Any help please?

This is what I got by @haci's code

```
> head(var.annovar.maf2[1:6,1:13])
   Tumor_Sample_Barcode Chromosome Start_Position End_Position Reference_Allele Tumor_Seq_Allele2 Hugo_Symbol Variant_Classification
1:                  new       chr1       89616151     89616151                G                 T        GBP7      Missense_Mutation
2:                  new      chr16       20556547     20556547                C                 G      ACSM2B      Missense_Mutation
3:                  new      chr17       66992092     66992092                C                 A       ABCA9      Missense_Mutation
4:                  new      chr20       13251339     13251339                C                 T        ISM1                 Silent
5:                  new       chr1          56987        56987               CA                 -           .                   <NA>
6:                  new       chr1        4133415      4133415               AG                 -           .                   <NA>
             tx   exon txChange aaChange Variant_Type
1:    NM_207398  exon6  c.C733A  p.L245I          SNP
2: NM_001105069 exon10 c.G1213C  p.G405R          SNP
3:    NM_080283 exon26 c.G3499T p.G1167C          SNP
4:    NM_080826  exon2  c.C327T  p.N109N          SNP
5:         <NA>   <NA>     <NA>     <NA>          DEL
6:         <NA>   <NA>     <NA>     <NA>          DEL
>

```

# Answer
> 1 votes
I am pretty sure that there is a reason that the first nucleotide of the `ALT` is the same as that of `REF` in `VCF` files, however I was not able to find a reference for this.

If you really need to do this for some reason, you will need to compare the alt and ref entries base by base starting from base 1 of each and until a mismatch. Here is some R code doing that:

```
library(data.table)
my_csv <- fread("example_vcf.txt")
print(my_csv)

# loop over each row
for(i in 1:nrow(my_csv)) {

  # make "strings" from the ref and alt columns, make "arrays" of letters out of these strings
  ref <- as.character(my_csv[i, "REF"])
  ref_array <- unlist(strsplit(as.character(my_csv[i, "REF"]), split = ""))

  alt <- as.character(my_csv[i, "ALT"])
  alt_array <- unlist(strsplit(as.character(my_csv[i, "ALT"]), split = ""))

  # which of the ref alt-pair shorter and how long
  shortest <- c(nchar(ref), nchar(alt))[which.min(c(nchar(ref), nchar(alt)))]

  # compare ref and alt base by base, remove from alt if same nucleotide present in both, stop at first mismatch
  index <- 1
  new_alt <- alt_array
  while(index <= shortest) {
    if (ref_array[index] == alt_array[index]) {
      new_alt <- new_alt[-1]
    } else {
      break()
    }
    index <- index + 1
  }

  my_csv[i,"ALT"] <- paste(new_alt, collapse = "")

  # use "-" if alt is empty
  if (my_csv[i,"ALT"] == "") {
    my_csv[i,"ALT"] <- "-"
  }
}

print(my_csv)

```

I have added some test cases at the end of the sample piece in your sample and the code is adding another column so that the changes made to the `ALT` are easier to see:

```
> print(my_csv)
    #CHROM       POS ID    REF            ALT       NEW_ALT
 1:   chr4   5708160  .      G         GATAAC         ATAAC
 2:   chr4  30280292  .      C        CACGTAT        ACGTAT
 3:   chr6 166253422  .      C        CCTTTCT        CTTTCT
 4:  chr11  21931340  .      A      AACACACGC      ACACACGC
 5:  chr12 118404422  .      G        GCACACA        CACACA
 6:  chr15  22672047  .      T TAAGTAGATTTAAG AAGTAGATTTAAG
 7:  chr16  32656293  . ATGATT              A             -
 8: chrXXX       XXX  .   TCTC              T             -
 9: chrXXX       XXX  .   TCTC           TCTC             -
10: chrXXX       XXX  .   TCTC             TA             A
11: chrXXX       XXX  .   TCTC            TAT            AT

```

---
Tags: vcf, text-processing
---


# Question
Title: which NCBI tool is optimized to identify a species from a DNA fragment?
“The Iceman” was a man who lived 5300 years ago and whose body was recovered from an Alpine glacier in 1991. Some fungal material was recovered from his clothing and sequenced.

Ice man : found as a frozen corpse in the Val Senales Glacier in Italy.

Thought to be carrying *Fomitopsis betulina* (formerly known as Piptoporus betulinus) commonly known as birch polypore, birch bracket, or razor strop) is a common Basidiomycota brown rot macrofungus growing on decaying birch wood.

Other species known : *Fomes fomentarius* (commonly known as the tinder fungus, false tinder fungus, hoof fungus,\] tinder conk, tinder polypore or ice man fungus)

Given two DNA fragments from NCBI genbank, try to find what species these DNA fragments belongs to or what species these fragments are most closely related to.

```
>Z54169.1 Iceman fungal clone 'Sim C-a1' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAAACTTGGGCCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTTCCTTCTGGAGAACCTCATGCCCTTCACTGGGTGTGTTGGGGAA
CCAGGACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAAT
AATAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCG
GGGGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCAT
TTGCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAG
TCTTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGA
GAAATCAAAGTTTTTGG

```

```
>Z54155.1 Iceman fungal clone 'Sim C-a2' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAACCTTGGGTCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTCCTTCTGGAGAACCTCATGTCTTTTATTAGGCGTGTTGGGGACC
AGAACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAATAA
TAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCGGG
GGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCATTT
GCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAGTC
TTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGAGA
AATCAAAGTTTTTGG

```

Now I know I can just google what species of mushrooms (which is where the above background info is from). But say that I was unable to and had to just go off the DNA fragments in order to find the species.

I've tried `BLASTn`, but I'm not sure if it is the right tool or if I need to properly calibrate the parameters to use it.

so far, I've tried:

Under “Choose Search Set:
Database : Non-redundant protein sequences (nr)
Organisms excluded: Animalia (taxid:33208), Plantae, Prokaryota  (taxid:2), Protozoa
Organisms included: fungi  (taxid:4751)
Included: Models (XM/XP),   Uncultured/environmental sample sequences
Under Program Selection:  Optimized for Highly similar sequences (megablast)

Algorithm Parameters (General Parameters, Scoring Parameters, Filters and Masking,
PS/PHI?DELTA BLAST)

are all kept default.
```

# Answer
My conclusion from the tree is the fungi appears to be genuine and not an extant contaminant (which is very easily done). It is not Fomitopsis betulina because this is division Basidiomycota. The fungi here is division Ascomycota and these are very distinct from Basidiomycota.

There have been many examples of contamination in ancient DNA work, but this doesn't appear to be one of them. The rationale is given below.

The first fragment was not interesting and described 100% identity to ascomycetes. The tree above is a simple distance tree from the second fragment of 18S which shows it is an outgroup to the ascomycetes (green terminal taxa). The ancient fungi contains 3 (I think) unique nucleotide residues that have not been sequenced in any extant ascomycetes to date and sequences are unique to any known species. The location of the 3 unique nucleotides is tightly clustered within the locus with 100% identities either side. This is significant because it demonstrates it is not a PCR mosaic artefact (yep its been done in ancient DNA work and subsequently proven artifactual). In the latter PCR mosaic artefact I'm referring to (which I wouldn't cite), was a fungi/bacterial mosaic from memory.

Given the homology of 18S in the extant ascomycetes I would say this appears to be a new genus, however you would need to be sure this was not sequencing error and this is very old DNA, which is subject to degradation even when frozen. The implication of the tree is quite unusual, it implies widespread extinction of this group of fungi in recent history (5000 years is slight in evolutionary terms) BUT it could be a sampling bias, i.e. no one has sequenced those fungi and NCBI is not representative.

I assume that you have performed phenotypic/morphological analysis on the spore to believe it is Basidiomycota.

The fragment in yellow is the query and the sequence next to it was the sequence already deposited on Genbank of this ancient isolate. You would need to look up the biology and ecology of ascomycete to understand its significance, but it is in a lineage of its own, so maybe difficult to interpret by comparison to extant species. You could concatenate the two alignments based on each fragment, but the problem with 18S is that is a hard to model and you'd need to perform secondary RNA structure modelling to assess whether these mutations would be viable.

In conclusion, it definitely looks genuine but it is not the species you think it is, in fact not even close.
> 1 votes

---
Tags: blast, phylogeny, phylogenetics, genbank
---


# Question
Title: What software can I use to develop a population haplotype network?
I looking for software to conduct a haplotype network using cytochrome oxidase (CO1) mitochondrial sequence data to study population structure in Anopheles mosquitoes.

# Answer
The software you are looking for is Network 5.0.0.3

Phylogenetic Network Software . Network © Copyright Fluxus Technology Ltd 1999–2018. Retrieved from https://www.fluxus-engineering.com/sharenet.htm

Whilst its old code it is still regularly used in pop gen for example there guys here recently used it. It will do the job, although personally I think Google Fusion tables are very cool, but that is complicated to use properly. R has some cool stuff in this context BTW, but in context I'd keep it simple (i.e. Network 5)
> 1 votes

---
Tags: phylogeny, phylogenetics
---


title: about GO terms's name

text:
# Question
Title: about GO terms's name
Did anyone know whether the GO terms can include more detail information? Like I can get the DotPlot of the GO terms as below:  The problem is that some of the genes, in the GO terms, is more about sinoatrial-node (SAN). So whether there is any method that we can do to dig deeper into the GO terms?

# Answer
> 1 votes
The various Gene Ontologies are arranged as Direct Acyclical Graphs. This means that each term can have children and parents. When a gene is annotated with a GO term, it also inherits all of that term's parents. Therefore, when you run a GO enrichment analysis, it is very common that you find the broader terms enriched since those are the ones shared by most genes.

For example, the Biological Process (BP) GO term GO:0003163 sinoatrial node development has the following tree:

```
￼ GO:0008150 biological_process
     ￼ GO:0048856 anatomical structure development
     ￼ GO:0032501 multicellular organismal process
         ￼ GO:0032502 developmental process
         ￼ GO:0007275 multicellular organism development
             ￼ GO:0048731 system development
             ￼ GO:0009888 tissue development
                 ￼ GO:0048513 animal organ development
                 ￼ GO:0072359 circulatory system development
                 ￼ GO:0060537 muscle tissue development
                     ￼ GO:0048738 cardiac muscle tissue development
                     ￼ GO:0007507 heart development
                     ￼ GO:0014706 striated muscle tissue development
                         ￼ GO:0003161 cardiac conduction system development
                             ￼ GO:0003163 sinoatrial node development
                                 ￼ GO:0060921 sinoatrial node cell differentiation

```

This means that any gene annotated with `GO:0003163` will also be annotated with every one of the terms above it in the tree, all the way to the top of the graph. All genes with at least one BP GO term will also be annotated with GO:0008150, the root of the tree.

As you can see, your "muscle tissue development (GO:0060537)" is also there and is a parent of "sinoatrial node development". However, if you just ran an enrichment analysis (I ave no idea what you did since you haven't actually told us anything), then what will have happened is that there are many genes annotated with various different child terms of GO:0060537, but since all child terms will inherit the parent annotations, what you see is the parent and not the specific term you are looking for.

So yes, there are far more specific terms that what you show, but whether or not they appear in your results will depend on what you are doing.

# Answer
> 0 votes
This ends up not being a question of what method to use, but rather how detailed the GO term database you're using is. GO terms are maintained by the Gene Ontology Consortium and are relatively broad for the most part, with a few more specific molecular signaling pathways included. If you search for the ontologies for "atrial" you'll find a number of relevant GO terms. The question mostly becomes whether those are sufficient for your needs. If they are, make sure you include "biological process" in your GO analysis. If not, you might have to make your own group of relevant genes.

---
Tags: rna-seq, go-enrichment
---


# Question
Title: Is there any tool for fuzzy sequence matching?
*Note: this question was also asked on biostars*

are you aware of any tool that is able to perform error-tolerant pattern-matching search on protein FASTA files?

For example, I want to know, which proteins in my fasta file to match`ADNG..C.G` regexp (which represents `ADNGCG` pattern). However, I want to be tolerant to matching errors, meaning that I'm good with any protein that differs in any 1 letter in the motif: `ADNG..C.D`, `ADNG..D.G`, `MDNG..C.G`, etc are all good. Running all possible variants through `grep` is possible, but exponentially long for longer patterns (I usually have `15-20` letters in pattern, and scan `~10e6-7` sequences).

I am aware of the tool `agrep` (docs. As far as I know, it's not supported anymore, and also does not distinguish letters in which it's allowed to make error. Also, it does not support long enough patterns (with more than 9 errors -- and yes, I tried to recompile it with myself from here and could not get it work). Also, it's not designed for proteins specifically.

I am also aware of ScanProsite and Protein Pattern Find services. The first one is not error-prone, the second one is seems to be a web interface for `grep`.

# Answer
Python solution:
e\<=2 will allow up to 2 insertions, 2 deletions, or 2 SNPs in any 2 combinations

```
import regex
seq_list=["GGGAGGG", "GGGAAGGG", "GGGAAAGGG", "TTTATTT", "TTTAATTT", "TTTAAATTT", "GGGG", "GGG"]
for sequence in seq_list:
    if regex.search("(GGGGGG|TTTTTT){e<=2}", sequence):
        print(sequence)

```

GGGAGGG
GGGAAGGG
TTTATTT
TTTAATTT
GGGG

Can also include a wildcard

```
seq_list=["GGGAAAAAGGG", "GGGAAAAATGG", "GGGAAAAATTG", "GGGAAAAATTT"]
for sequence in seq_list:
    if regex.search("(GGG.*GGG|TTTTTT){e<=2}", sequence):
        print(sequence)

```

GGGAAAAAGGG
GGGAAAAATGG
GGGAAAAATTG
> 2 votes

---
Tags: proteins, python, blast, software-recommendation
---


# Question
Title: move files of a different group to separate directory
> EDITED

I have list of my samples in first column and corresponding file name for my .vcf files in second column of a file name `clin_name.txt` like below. For example name of one vcf is `LP6008031-DNA_A04_vs_LP6008032-DNA_D03.snp.pass.vcf`

```
OC/AH/120/E/T3  LP6005409-DNA_F01
OC/AH/126/E/T1  LP2000325-DNA_A01
OC/AH/136/E/T3  LP6005409-DNA_E02

```

This is the link of my file https://www.dropbox.com/s/mvzf47feiqe7x69/clin\_name.txt.txt?dl=0

In the second file name `pre_post.list` I have information of pre or post treatment for each sample like this

```
pre-treated post-treated
OC/AH/120   OC/AH/065
OC/AH/126   OC/AH/069
OC/AH/136   OC/AH/074
OC/AH/160   OC/AH/084

```

This is link of this https://www.dropbox.com/s/nhv3odjkdjdp7hz/pre\_post.list.txt?dl=0

For moving different .vcf files based on the pre or post treatment to separate folders I typed

```
cat pre_post.list | awk '{print $2}' | sed 's@.*/@@'  > post_list.txt
cat pre_post.list | awk '{print $1}' | sed 's@.*/@@'  > pre_list.txt

```

Then I have

```
pre-treated
120
126
136
160
182

post-treated
065
069
074
084

```

I then typed `./test.sh` in which I wrote

```
for item in  $(cat  pre_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/
    else echo post_list $item;  fi;  done;

```

That returns

```
cp: cannot stat `*OC/AH/069/S/T3*': No such file or directory
[fi1d18@cyan01 snp]$

```

Or I run

```
for item in  $(cat  post_list.txt | grep -Eo '[0-9]{1,4}');
    do if grep -q $item clin_name.txt ;
    then out=$(grep $item clin_name.txt | awk  '{print $1}');
         cp *$out* /temp/hgig/fi1d18/TRG45/snp/pre/
    else cp *$out* /temp/hgig/fi1d18/TRG45/snp/post/;  fi;  done;

```

That returns

```
cp: cannot stat `OC/ST/026/S/T1*': No such file or directory
[fi1d18@cyan01 snp]$

```

I really don't know how to solve this :(

# Answer
You need to bring together the information in the two tables, I prefer `R` for such actions. The common columns between the two tables do not quite contain the same information hence additional data cleaning was required. Moreover the table with treatment information had to be changed to long format with `melt`. And I have filtered the `merged_data` so as to conation only "pretreatment" info but you can add a if/else statement to the for loop to control how/where you will move your files.

```
library(data.table)

clin_name <- fread("clin_name.txt", header = TRUE)
clin_name$sample_id <- sub("/E/T.", "", clin_name$sample_id)

pre_post <- fread("pre_post.list", header = TRUE)
print(pre_post)
pre_post <- melt(pre_post,
     measure.vars = c("pre-treated", "post-treated"),
     variable.name = "treatment",
     value.name = "sample_id")

> print(pre_post)
      treatment sample_id
1:  pre-treated OC/AH/120
2:  pre-treated OC/AH/126
3:  pre-treated OC/AH/136
4:  pre-treated OC/AH/160
5: post-treated OC/AH/065
6: post-treated OC/AH/069
7: post-treated OC/AH/074
8: post-treated OC/AH/084

merged_data <- merge(clin_name, pre_post)
print(merged_data)

pre <- merged_data[merged_data$treatment == "pre-treated", "file_name"]

> print(merged_data)
   sample_id         file_name   treatment
1: OC/AH/120 LP6005409-DNA_F01 pre-treated
2: OC/AH/126 LP2000325-DNA_A01 pre-treated
3: OC/AH/136 LP6005409-DNA_E02 pre-treated

# invoking a terminal command to move each file to a directory called "pre"
for (i in pre$file_name) {
  command <- paste0("mv ", i, " pre/")
  system(command)
}

```
> 1 votes

---
Tags: shell
---


# Question
Title: What is the state of the art for GWAS in terms of the statistical algorithm for either Case/control and Quantitative traits?
*This question was also asked on Biostars*

I'm trying to understand what is the best algorithm for GWAS nowadays. I know we have many tools available like Plink and Hail, but currently, what is the best algorithm if I won't use any them? Let's say, write down a script in R or Python from scratch. Which statistical algorithm should I use? Is it Linear mixed models (LMMs)? I'm confused as we can have binary phenotypes (case/control) or quantitative phenotypes. LMM seems to address quantitative ones, but can it be used for case/control as well? Actually, what is the state of the art for both/each of them? Pair-reviewed papers as references will be appreciated. Thanks!

# Answer
> 2 votes
I've got the answer to my question at https://www.biostars.org/p/410466/ by chrchang522

> "The main regression executed by Plink was introduced by EIGENSTRAT in ~2006; see https://www.nature.com/articles/ng1847 . This is actually straightforward to write in R/Python from scratch; the harder part is optimizing the implementation for large datasets.
>
> The Firth regression added to Plink 2.0 to improve handling of rare variants and imbalanced binary phenotypes was motivated by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4049324/ .
>
> Mixed linear models provide better statistical power when you have lots of close relatives in your dataset, but are much trickier to solve; actually, this is still a significant research area. Two tools covering parts of the current state-of-the-art are SAIGE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6119127/ ; handles imbalanced binary phenotypes, but relatively slow) and fastGWA (https://www.nature.com/articles/s41588-019-0530-8 ; great speed, but doesn't support dosage data yet and uses a misspecified model for binary phenotypes)."

---
Tags: statistics, gwas
---


# Question
Title: Which of these enzymes is a better target for virtual screening?
I am a MSc student in Applied Biochemistry. For my Masters Thesis, I will be doing some research about finding new potent inhibitors for a specific enzyme of my choice using molecular docking (virtual screening).

I have selected 5 enzymes that have all been approved by my advisor and now the choice is for me to make. I am very interested in neurological diseases as you can see:

Kynurenine 3-monooxygenase (KMO)

Glycogen synthase kinase-3 (GSK3)

Dopa decarboxylase (DDC)

HMG CoA reductase

Phospholipase A2 (PLA2)

Knowing that I plan to publish my work. I would be really really grateful, if you guys could help pick :)

# Answer
The ultimate choice stems from your individual interests and not a scientific reason... but having said that some pointer of issues may be handy.

I am assuming that you want to make an inhibitor that messes with the activity and not a ligand that binds somewhere and you then plan to add a moiety to make it a protac to get it degraded or to make it fluorescent to mark it. If not, ignore any talk of activity.

## Crystal

* is a crystal structure available?
* does the structure have the cofactor if any?
* is a structure available (itself or homologue) with the native metabolite bound?
* how is the active site? In the case of many cytochrome P450 enzyme the entranceway is really deep and is guarded by a flappy gate, both of which confer a great deal of specificity regardless of the residues in the active site, which are often just water coordinating. Your first enzyme, is that a NADPH/FAD or a cytochrome enzyme?
* What is the mechanism and is there a classic type of inhibitor chemistry? I.e. catalytic triad with a cysteine nucleophile is a nice target.(1)

## Methodology

* Wet lab: A problem with computational chemistry/biochemistry is that someone needs to make the protein and physically test your compounds. Are these easy to express? This may sound like not your problem, but a lot of issues can come from this latter stage.
* is there a close homologue you want to avoid targetting? If so, you'll have to do twice the docking (not a big deal and potentially more interesting). (2)
* is the method you plan to use okay with covalently bound ligands? Dopa decarboxylase sounds like a PLP enzyme. If it is, making an inhibitor using classic functional groups that nucleophilically attack the pyridine ring is very straighforward, but not all docking methods are okay with covalent ligands.
* Are any hits likely to follow Lipinski's rule of five? If you have a very charged active site, your hits will be metabolite-like (polar) and not drug-like (horrendously insoluble).

A minor thing is to make sure that nobody else is working on it or you may be left behind. Industry is all cloak and dagger, but the academic SGC has a list of \[target enabling packages)\[https://www.thesgc.org/tep\] and their future ones (nominations).

## Physiology

This was most likely in your lit search. But I am double checking.

* You have chosen these enzymes because their unchecked activity is disease-causing, right?
* And they do not have any other roles in other tissues?
* Abolition of activity does what you want?

## Other hits

Lastly, but not least. See which of these have many paper about barely active hits: this should mean it is hard. For example, bacterial alanine racemase has been a holy grail for a while, but progress is slow as it has been described as a "ship in a bottle" due to its long and unforgiving entranceway.

## Addenda

1. A catalytic triad is acid-base-nucleophile: if you have an electrophile with a leaving group (including internal ones, such as epoxydes) you can virtually screen more complex small molecules that resemble your ligand, but with the difference that there is a leaving group, resulting in an irreversible inhibition. Whereas competitive inhibitors require nanomolar affinity to outcompete the natural ligand, an irreversible inhibition can even slightly worse and will work.
2. The problem of homologues.
    * Paralogues (homologous genes in the same organism): you want to target protein 1, but you do not want to target protein 2. This is a common problem with kinases. You therefore need to find a hit that works better in one than in the other, which with a big of rational inference is easy.
    * Homologues (homologous genes in different organism). You want to kill species A (say a fly), but not B (say a human). Same thing.
3. Twice the docking. It is more of a challenge, but makes for a better spin on the story. There is an explosion of people doing virtual screens that do not go anyway afterwards, so it is getting harder and harder to publish studies without extensive empirical backing. Hence the need for an interesting angle.
> 1 votes

---
Tags: docking, molecular-genetics
---


# Question
Title: Which taxonomy data type to use for abundance graphs and statistical analyses?
At our institution we have had some debate about this topic. So let me give some background:

Qiime2 produces a OTU table which shows the number of sequences that matched to a specific OTU. It also produces an Taxonomy table which links all the OTUs to their corresponding taxonomic identity.

So what we were discussing was whether to use the number of sequence hits per OTU for alpha and beta diversity, as well as abundance graphs. Or should we use the Taxonomy table and count the number of occurrences for each taxa.

An example of this would be something like the following:

Here are two OTUs:

`D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Chloroplast;D_4__Virgulinella fragilis;D_5__Virgulinella fragilis;D_6__Virgulinella fragilis` and `D_0__Bacteria;D_1__Cyanobacteria;D_2__Oxyphotobacteria;D_3__Nostocales;D_4__Oscillatoriaceae;D_5__Oscillatoria PCC-10802;Ambiguous_taxa`

If we used the first approach then there would be, for example, 3000 sequence matches for the first OTU and 1500 for the second OTU. However, they represent two different species from the same Class. So abundance wise, the first OTU could be represented twice as much as the second. With the second approach they would share the same abundance percentage.

What is the correct approach for these types of analyses according to the bioinformatics community?

# Answer
> 0 votes
The fundamental answer to this question is the role of allele frequency, i.e. population genetics, within taxonomy. The answer of course is that population genetics and population dynamics should be a integral part of taxonomic designation PROVIDING it is not subject to unreported sample bias. Thus it has to be option 1.

However, your definitions I'm not so sure about because the counts per OTU using your literal definition means every SNP within a species. I don't believe that is what you mean, but that is what you have said. Obviously at a taxonomic focus of a single SNP this isn't meaningful because there would be as many OTUs as there are sequences and the species definition would be meaningless.

Both approaches suffer from lacking a clear definition to make a species assignment and this is controversial. If you are using a tree to assign species membership, e.g. 16S at 5% identity, then that would generally be acceptable, but some would question it.

You can of course use deep learning algorithms to do this, but this is only as accurate your training and at present (check Bioinformatics journal) the accuracy is at best 90%, which is very low by comparison to phylogeny.

Thus I think it is a non-question, because allelic diversity is extremely important in vast swathes of biological analysis. Taxonomy is not a disapline in its own right, but a foundational step biologists use to understand the ecotope, ecoregion etc... The only justification for pursing approach 2 is if there is persistent concern in the bias of the sampling process, to the point where the alleleic (species) frequency becomes meaningless.

---
Tags: taxonomy
---


# Question
Title: Interpreting quantitive outputs from maximum likelihood phylogenetic trees
I ran a calculation in RAxML to determine the majority consensus phylogeny of a maximum likelihood bootstrap (How to show bootstrap values on a phylogenetic tree constructed with RAxML), and I got three output files:

1. ***RAxML\_bipartitions.output\_bootstrap.tre***
2. ***RAxML\_bipartitionsBranchLabels.output\_bootstrap.tre***
3. ***RAxML\_info.output\_bootstrap.tre***

What is the difference between the

* ***RAxML\_bipartitions.output\_bootstrap.tre*** file and
* ***RAxML\_bipartitionsBranchLabels.output\_bootstrap.tre***?

# Answer
> 1 votes
In summary,

```
RAxML_bipartitions.output_bootstrap.tre

```

Is the only file of interest. The reason this is true in this context is really complicated and you have to understand the statistics of likelihood and how they are interpreted within phylogeny to understand why. This file is simply the final output of a non-parametric bootstrap analysis performed by maximum likelihood.

**What on earth is a non-parametric boostrap?**

A non-parametric bootstrap is resampling each alignment position with replacement. Thus if we have alignment positions 1,2,3,4,5 A bootstrap resample for 2 replicates might be,

Replicate 1

```
1,1,3,5,2

```

Replicate 2

```
4,2,5,2,1

```

The ML algorithm will make trees of replicates 1 and 2 and find the consensus between them. If you think about it in any other context a bootstrap replicate is pretty meaningless because it no longer reflects the true biological sequence. Thus information on how the consensus was derived, are not really of interest to us providing we are confident this has been done correctly, viz. RAxML\_bipartitionsBranchLabels.output\_bootstrap.tre and RAxML\_bipartitionsBranchLabels.output\_bootstrap.tre

**So why is this output of limited use?**

There are situations to some investigators this information is useful, but assess the robustness of a tree topology its not needed. The only thing we want is a phylogram (bestTree) with the bootstrap values superimposed on them. We really don't need complicated stuff such as the tree to be represented for example as a polytomy (non-bifurcating tree) because we can just read the bootstraps to make that deduction (values \>\> 75%). In addition, there is not perfect consensus what boostrap value constitutes robustness, but generally most agree \>80% is robust.

---

**What output files have useful information in them?**

The information that is important are the files associated with "bestTree", that was the single maximum likelihood tree performed on the intact native sequence. The "info" file for this contains 3 really important parameters:

* -lnL ... very important!!
* Gamma distribution parameter "alpha",
* PINVAR, proportion of invariant sites,

**-lnL** is the highest log-likelihood (probability) of the phylogeny. It is usually a very small number for which where is an enormous amount of theory over it.

**Alpha parameter of the gamma distribution** this is the shape parameter of the mutation rate, if it is very low (\<1) the distribution of mutations across the alignment is very tight clustered and approximates to a negative binomial distribution. This means some sites don't mutate at all and a small number of sites mutate alot. If it is very large \>200 (which is never observed) it approximates to the Poisson distribution, meaning the mutation distribution is randomised across the alignment.

**PINVAR** this is a straight percentage/frequency and simply means the sites that don't mutate.

**How are they calculated?**

PINVAR and alpha are not emperically calculated, i.e. if you look at an alignment and say 'no mutations at that position', PINVAR would of course agree but may consider other invariant depending on the phylogeny. These parameters are calculated by maximum likelihood and you can begin to see why the calculation takes so long ... alpha and PINVAR affect the tree topology (which affects -lnL), but the topology affects alpha and PINVAR. Thus, is a multidimensional search of tree and parameter space.

**So what stuff do I report in my Results?**

Anyway reporting -lnL is good technique and shows the reader you've done maximum likelihood, citing PINVAR and alpha from gamma distribution helps ('Methods' parameters were calculated reiteratively under maximum likelihood). This is only useful for bestTree. The -lnL, PINVAR and gamma's alpha are also calculated for every single bootstrap replicate, but these values are of limited to use, because we have resampled the data, only the consensus tree counts... Obviously presenting the bootstrapped phylogram is extremely important.

Welcome to the technical world of phylogeny!

The amino acid matrix you used BTW .. LG is in vogue right now.

**How do I do it?**

When I do this stuff its via Biopython and ETE3, I capture the values within the pipeline and don't examine the output files of RAxML because I generate my own.

---
Tags: phylogenetics, phylogeny, software-usage
---


# Question
Title: how to interpret of "pvclust" dendrogram and finding height for cutting dendrogram?
I study on RNA-seq expression dataset about one cancer in TCGA. I downloaded FPKM dataset and removed batch effect by ComBat() function. Now, I used `pvclust` for clustering 72 samples and the result of that can access in this link.my code for running pvclust is as below:

```
fit_LMS <-pvclust(datForpvclust,method.hclust = "average",
                          method.dist = "euclidean",nboot=10000, parallel=TRUE )
pvrect(fit_LMS, alpha=0.95)

```

As you see, by alpha=0.95 just 3 clusters are significant although in right cluster a few clusters have been included. really I can't analysis this result. So, I have 2 questions: 1- What is the interpretation of `pvclust` dendrogram? Does my dataset has meaning full clusters?

2- I am interested in the height of tree cut in the dendrogram, which height is better for this dataset based on pvclust result? H=105 or H=110 or another height?

I appreciate it if anybody shares his/her comment with me.

# Answer
> 1 votes
* You do not need to give links to your images, you can embed images in your questions using the "image" icon.
* FPKM (and RPKM) might not be the best RNA-seq normalization method, see this and this.

`pvclust` providing p-values for the clusters computed by `hclust`. To do so, bootstrap resampling is used. From their documentation:

> Values on the edges of the clustering are p-values (%). Red values are AU p-values, and green values are BP values. Clusters with AU larger than 95% are highlighted by rectangles, which are strongly supported by data.

Basically you are getting a clustering pattern and a statistic associated with each of the clusters giving you an idea about the "statistical significance" of the clusters you have obtained (i.e. clusters *most probably* not being a result of strong batch effetcs, ...).

To be able to comment on the resulting clusters making sense or not, you will need to consult somebody with *domain knowledge*, it is not a number or two that would determine the quality of the output but also how good is the fit between your result and the (expected) biology. For example in your case with TCGA data:

* Do the samples cluster by batch?
* Do tumor and adjacent normal tissue samples aggregate into different clusters?
* Do the samples in the same cluster share the same mutation/methylation pattern/disease subgroup?

---
Tags: r, rna-seq, clustering
---


# Question
Title: Tree cut height for Consensus data input not working
I have tried different combination of tree cut heights but im getting the exprSize as one sample as full the other one as empty .

```
# Choose the "base" cut height for the female data set
baseHeight = 135
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(135, 145*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

```

My figure

When i see the cutHeights and exprSize i do see its higher than 120 but still i get only one set of sample which is 55 whereas the my other set which is

```
exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55 47

$structureOK
[1] TRUE

```

This is after tree cut

```
exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55  0

$structureOK
[1] TRUE

> cutHeights
[1] 135.0000 123.9091

```

Im not able to figure out what is wrong

# Answer
Your data don't have obvious outliers so you don't have to worry about removing them. The problem is that both sets show prominent sample clusters. I would try to figure out what drives those clusters (is it technical or biological), and whether the variation should be removed. See the comments in WGCNA FAQ under point 5 (My data are heterogeneous. Can I still use WGCNA?). This looks like data from GEO, so look into the sample annotation for these data sets and see if any of the recorded sample characteristics align with the clusters you see in the sample trees.
> 1 votes

---
Tags: wgcna
---


# Question
Title: Are crystal structure of Ligands used in volunteer computational docking like FightAIDS@home?
In large Volunteer computing for molecular docking projects such as FightAIDS@home, do they use the original crystal structure of ligands in docking or unconfirmed structures of ligand?

# Answer
I don't know that project's details, so cannot say officially, but undoubtedly it will be

* crystal structures of protein or models of point mutants based on structures of protein.
* Computer structure of the ligands based on force-field calculation

## Protein structures

Firstly, HIV virus has been extensively studied by structural biologists, so most components have structures.

Secondly, docking on models is generally done only for threaded models based on closely related templates, say a human model of a mouse structure or a recently diverged pair (i.e. conserved proteins that split with the genome double duplication in bony fish). Whereas a given virus protein will likely have poor homology to a known structure and not be a good candidate.

A corner case is core point mutants... If you have a mutation in a protein it is likely going to be on the surface and not core: a model of a point mutant is likely very faithful to an empirical version —and generally you'd want to design drugs to a region that is not variable anyway as you'd want it to work universally. If you have a change in the active site to confer resistance say (i.e. in the core), there is generally a fitness loss to the organism (say lowered catalysis or thermostability) which needs to be compensated with other mutations (epistasis). As this is rarer and very interesting scientifically and medically, these resistant variants will most likely have been crystallised as the conditions for crystallization will most likely be the same.

## Small molecule structures

Force field calculations for small molecules are most often valid (there computer Vs empirical RMSD is minute). The partial charges may be off, but that's not something crystals give anyway: In fact a ligand used in docking is (generally/ought to be) energy minimised beforehand. For docking getting charges wrong is the killer.

Crystallization of small molecules is used as validation of a synthesis route, so some may have structures, but they will not be used in a pipeline as a SMILES string is all you need.
> 1 votes

---
Tags: docking
---