# Question
Title: Algorithmic investigation of viral proteins for vaccine design
I wish to investigate how a virus could be attenuated through targeted mutagenesis by a priori prediction of attenuating amino acid mutations.

In particular, I am seeking bioinformatic solutions where the protein sequence of an exisiting live vaccine strais has been compared against it wild-type counterparts to identify key attenuating mutations and the key attenuating mutation(s) can be deduced without wet-lab site directed mutagenesis.

I presume surface antigens, such as the spike protein in COVID-10 would be the proteins of interest, but I'm open to suggestions.

I am a computer scientist looking to investigate the mutational parameter space of viral protein sequence data. A key approach I am considering is to use a mutational correlation/covariance matrix and then remove convolution artefacts. I hope this will identify compensatory mutations that will maintain a local interaction for example within a fold, providing a starting point for identifying attenuating mutations.

This may have application for example to COVID-19 and the goal is to produce an efficacious vaccine.

# Answer
> 5 votes
There isn't a vaccine for any coronavirus, and your question is generally about targeted attentuation, which is a complex area.

The basic building blocks for any vaccine development is virological understanding of the proteins involved in pathogenesis. I will focus on covid-19 as an example here.

The majority of bioinformatics work is based around the spike protein as is wetlab vaccine development. This gene is immediately after the 20 kb of ORF1ab, including the RNA dependent RNA polymerase, its around 3.5 kb but it might be a bit bigger. The spike protein is a major B-cell (antibody) antigenic determinant and the cellular receptor sites are located within it. Your bioinformatics could ultimately focus on targeting these specific sites, in particular using antigenic (antibody) prediction algorithms.

First you need some basic biology of spike and hydrophobicity prediction is a good starting point such as TMHMM here and consider the SARS spike protein structure for example here. Once you've a handle on those structures antigenicity scanning would be a next phase analysis and possibly homology modelling. This is because COVID-19 and SARS are closely related, but there are key differences. You finally consider the full range of amino acid mutations in current epidemic via NCBI.

At this point you would target a specific vaccine design strategy and have a clearer understanding on the approach you wish to take. Good luck!

---
Tags: sequence-alignment, phylogenetics, phylogeny, covid-19, sars-cov-2
---


title: Practical use of RNA structure

text:
# Question
Title: Practical use of RNA structure
I have spent some time studying RNA structures. While there is a range of interesting phenomena and functions that certainly deserve understanding from the scientific point of view, I have never encountered a practical use of this knowledge (in terms of disease diagnostics, medications, etc.) Have you? If a practical use of RNA structures is still in the future, what could that be?

# Answer
> 3 votes
Please look up flavivirus 'double loops' as you described them previously (post for "Coronavirus RNA') and associated RNA secondary structure anomalies for dengue virus and associated vaccine (Butantan) and the yellow fever virus and its vaccine (17D). If you are aware of "double loops", we must be aware of the association of RNA secondardy structure and attenuated strains, i.e. live vaccines. The point is that modification of the RNA secondary structure is a key strategy in flavivirus attenuation.

---
Tags: phylogenetics, phylogeny, covid-19, rna-splicing, rna-structure
---


# Question
Title: How to seperate the genus into different group based on th ANI calculation?
I've calculated the ANI of the streptomyces genus, but I'm confused about how to separate them into a different group, could anybody kindly tell me how to do that? Thanks!!!

# Answer
> 1 votes
ANI = average nucleotide identity.

This will give you a 2x2 matrix. You reconstruct a phylogeny using clustering, the most important being neighbor-joining, UPGMA is very occassionally used.

The statistic is clustering algorithm and its done by pairwise reconstruction of sister groups. Once a sister-group is defined (usually on nearest-neighor both in neighbor-joining and UPGMA) the genetic divergence from the hypothetical node is established either using the mean (UPGMA) or additive criteria (neighbor-joining). Okay enough theory ...

You get Mega X and find out the import format for a 2x2 distance matric. Mega X however will calculate this for you. You then go to the build tree menu heading and click "neighbor-joining" (or "distance"). The tree will very quickly automatically appear.

---
Tags: phylogenetics, phylogeny
---


# Question
Title: What does the yellow mean in this image from Virus Pathogen Resource?
What does the yellow mean in this image of a SARS (coronavirus) protease from the Virus Pathogen resource database? It's from https://www.viprbrc.org/brc/

# Answer
Sulfur atoms are shown in yellow.

The molecular viewer that you use is JSMol (JMol ported to the web).

Atoms are colored by element: grey C, blue N, red O and yellow S.
If you wonder how other atoms would be colored, see JMol's Default element colors, by periodic table.
> 7 votes

---
Tags: proteins, protein-structure, covid-19, sars-cov-2, molecular-genetics
---


title: Is it possible to do this in R?

text:
# Question
Title: Is it possible to do this in R?
I have a big clinical data file with 96 columns like age, gender, BMI, etc

I want to see which of these clinical characteristics respond to chemotherapy. Response to chemotherapy divides patients to two groups as `yes` and `no`. Among these clinical characteristics, some are based on continuous data and some reflect \`categorical data. Is it possible to write a function to loop over these characteristics and test for relationship between response to chemotherapy? Something like Wilcox test for continuous values and chi square for categorical characteristics. By hand, it takes too much to do for all one by one !

This is my data

```
ID  response_to_chemo   BMI Chemo   Predictor   DJANGO  gender
AH/155  no  50  1   1   0   M
RS/022  no  67  1   1   0   M
RS/027  no  80  1   1   1   F
ST/023  no  65  1   1   0   M
SH/051  yes 47  1   1   0   M
AH/075  yes 90  0   1   0   M
RS/047  no  67  0   1   0   F
ST/029  no  61  1   1   0   F

```

For instance

```
myTable <- table(clin$gender, clin$response_to_chemo)

chisq.test(myTable)

```

Gives p-value for testing the relationship of gender (CATEGORICAL) in response to chemotherapy

And

```
t.test(clin$BMI ~ clin$response_to_chemo)

```

is for BMI

I meant a function calculating these p-values for clinical characteristics one by one

# Answer
> 6 votes
R supports logistic regression, which would seem to be the most efficient method for tackling this question. Assuming the "Chemo" variable is the type of chemo the code would be something like:

```
glm( (response_to_chemo == "yes") ~  BMI + Chemo  + Predictor +  DJANGO + gender,
        family="binomial")

```

EDIT: Corrected a typo (added a missing double quote)

# Answer
> 3 votes
The general approach would be to loop through every column starting at column 2. You can use numeric indexes to do that.

For each column, check its type. If it is a factor, use your `chisq.test` method. If it is of numeric type, use your `t.test` approach. Write the results to a list so you can parse them later.

# Answer
> 2 votes
You have category-based variables as well as variables based on continuous data. The best-option could be linear Mixed- effects model for your research purpose. Chi-square test and t-test may be used to cross-validate or re-examine the outcome of LME model.

---
Tags: r, database, statistics, software-recommendation, correlation
---


title: Coronavirus RNA structures?

text:
# Question
Title: Coronavirus RNA structures?
Is there anything known about the RNA structures of coronaviruses? More specifically - do they have any interesting known structures in the translatable region, like RRE of HIV or the double loops in flaviviruses?

**Update**
Here is a recent development on the side of structure prediction.

# Answer
> 3 votes
To add a more complete answer: the current coronavirus is closely related to the SARS virus that caused the outbreak in 2004, and on which much research has been done.

Here is a general review of the coronavirus epidemiology, life cycle etc.

I haven't found yet any materials about the RNA structures in the translatable region, however the structures in the 3' and 5' terminal regions are well studied and are known to be crucial for initiating viral replication and for its packaging.

This recent article addresses in details the RNA structures in 5' terminal of human coronaviruses, whereas this older article describes the structures both in 3' and 5' terminals.

Finally, the interesting algorithmic developments for the link that I added to the question are described here.

---
Tags: rna, covid-19, rna-splicing, sars-cov-2, rna-structure
---


# Question
Title: Are there any Sulfur Dioxide molecules in this structure?
Does the protein structure shown in the JSmol link have any Sulfur Dioxide molecules?

In the picture, the Sulfur element is yellow and Oxygen is red (see here)

https://www.viprbrc.org/brc/structureOperation.spg?accession=1P9S&decorator=corona&context=1583712525175

Link to mmcif file (I have the VS Code plugin installed):

https://files.rcsb.org/view/1P9S.cif

# Answer
The answer, as provided in the comments by gilleain is No. The commenter wanted to know why I asked this. The reason is I heard sulfur dioxide as silent speech when thinking about antibodies that’s why I asked.

I searched the file https://files.rcsb.org/view/1P9S.cif for the SO2 molecule and didn't find it. I got confirmation in the comments of the question.

Would anyone be interested in an app that takes a .cif file of a molecule and allows you to search the file for any type of element or molecule in it? Would that help anyone with their research?
> 2 votes

---
Tags: protein-structure, molecular-genetics
---


# Question
Title: How to get the intersection of peaks after peak calling using MACS2?
I have following 5 narrow peak files after peak calling.

```
K14_peaks.narrowPeak
K15_peaks.narrowPeak
K16_peaks.narrowPeak
K3_peaks.narrowPeak
K8_peaks.narrowPeak

```

I need to get the intersection of peaks into a new file (The peaks that represent all these files). I used following code but I got higher number of peaks not a lower number.

```
bedtools intersect -a K14_peaks.narrowPeak -b K15_peaks.narrowPeak K16_peaks.narrowPeak K3_peaks.narrowPeak K8_peaks.narrowPeak -wo > intersected.bed

```

Can someone help me? Thank you.

# Answer
> 1 votes
`bedtools intersect` with multiple files after `-b` performs pairwise intersections between `-a` and each file in `-b`. It does not intersect all files simultaneously. See here, from the Bedtools documentation. The green section only contains intervals from `-a` that also intersect at least one interval in any `-b` file.

What I think you're looking for is `bedtools multiinter` (documentation here is unfortunately blank). `bedtools multiinter -i input [input...]` will produce a BED file that covers every position listed in your input files, followed by how many and which of your input samples overlap that given position.

You can use this BED file and create something like an UpSet plot to visualize how many bases/peaks are shared between any combination of your samples.

---
Tags: ngs, bedtools, chip-seq, peak-calling, macs2
---


# Question
Title: processing multiple fastq files with cutadapt
I have DNA sample from 5 pools, having 25 fastq files each. I am running cutadapt to remove the primers using this command `cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG A7_S7_R1_001.fastq -o sb_A7_S7_R1_001.fastq --discard-untrimmed` However, this is too time consuming. I tried my hand at writing shell script

```
for file in /dir/*
do
  cmd cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

```

But it is not running . Can someone help. Thanks

# Answer
> 2 votes
Try this one

```
for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

```

Also you have omitted `--discard-untrimmed` from the original command. If you want that as well along with everything to run in parallel:

```
for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG --discard-untrimmed "$file" -o "`basename -s .fastq $file`_out.fastq" &
done

```

---
Tags: fastq, shell, cutadapt
---


# Question
Title: How to interpret Mendelian randomization results?
I did Mendelian randomization using this software: https://cran.r-project.org/web/packages/MendelianRandomization/vignettes/Vignette\_MR.pdf

```
library(MendelianRandomization)
f=read.table("246LDout272Biobank_Retina.txt", header=T)
> head(f)
          rs exposure.beta exposure.se outcome.beta outcome.se
1  rs1029830      0.723525  0.03026430  0.066715400  0.0359278
2  rs1029832      0.723785  0.03029603  0.064105600  0.0359021
3 rs11078374     -0.411789  0.04189295 -0.000376929  0.0406439
4 rs11078382      0.882549  0.14275799 -0.197074000  0.1247720
5  rs1124961     -0.333763  0.05589377 -0.075468600  0.0576012
6  rs1135237     -0.316831  0.05552530 -0.074086200  0.0573111

MRInputObject <- mr_input(bx = f$exposure.beta,bxse = f$exposure.se,by = f$outcome.beta,byse = f$outcome.se)
EggerObject <- mr_egger(MRInputObject,robust = FALSE,penalized = FALSE,correl = FALSE,distribution = "normal",alpha = 0.05)

MR-Egger method
(variants uncorrelated, random-effect model)

Number of Variants =  246

------------------------------------------------------------------
  Method     Estimate Std Error  95% CI       p-value
  MR-Egger    0.115     0.016  0.084, 0.146   3.28e-13
 (intercept) -0.010     0.008 -0.025, 0.006   0.226
------------------------------------------------------------------
Residual Standard Error :  1.106
Heterogeneity test statistic = 298.4508 on 244 degrees of freedom,     (p-value = 0.0099)
 I^2_GX statistic: 97.2%

```

Let's say that my exposure is called Retina and my outcome is called Biobank. Can someone please help me interpret these results in terms of 'horizontal pleiotropy' as mentioned here: https://www.ncbi.nlm.nih.gov/pubmed/29771313

# Answer
> 1 votes
It looks like you have linkage as opposed to linkage disequilibrium, therefore horizontal pleiotropy is an explanation of the calculation.

The calculation appears to be a GLMM around the beta distribution, the residual looks low (which is good). The heterogeneity statistic is significant, you would need to assess whether this is an index of heterozygosity, but it will relate to expectation of Mendellian inhertiance based around the Hardy Weinberg equilibrium.

---
Tags: r, phylogenetics, phylogeny, statistics, population-genetics
---


# Question
Title: What is the significance of the effect size on WOMAC scores for TKA?
https://www.medicaljournals.se/jrm/content/html/10.2340/16501977-2210

I understand the difference between statistical significance and effect size, but I am having trouble interpreting the data on Table III in the article linked above. I am evaluating the WOMAC scores. My first question is can effect size confidence intervals cross zero and still be valid? My second question is about the negative and positive values of the Hedges'g. The first study scores the WOMAC in a way where the higher the score, the better the outcome. The second two studies score the WOMAC in a way where the lower the score, the better the outcome. Does this mean that in the first study a positive Hedges'g indicate larger effect for the treatment group and a negative Hedges' g for the other two studies indicate larger effect for the treatment group?

# Answer
> 1 votes
Hedge's g measures the difference between two groups here experimental and control in terms of combined (weighted) standard deviation. A higher g implies greater difference and lower g means lesser difference between two groups. It depends on the study how the scoring of questions is proceeded with : the effect-size will be same. For example, answer to question item may be scored 1 to 10 depending on whether question item is written in positive sense (the apple is good for health) or negative sense (the apple is not good for health). Positive Hedge's g does not indicate larger effect for the treatment group. THE negative Hedges' g for the other two studies indicates positive effect for the treatment group if the difference has been computed by deducting experimental group mean from the control group mean (which is less than exp. Group Mean). Overall, I find that a knowledge of questionnaire scaling and scoring of scales and question items could be very helpful in general and in doing a meta-analysis effectively.

---
Tags: statistics
---


title: Pausing index calculation

text:
# Question
Title: Pausing index calculation
I have Pol\_II data for Ctrl vs Knockout condition with two replicates each. I want to calculate the pausing index for the Ctrl vs KO condition.

Here I have given tag directory as input to Homer.Now what I see is as such,for each biological replicate it calculates the pausing ratio.How can I calculate as a whole like 2 Ctrl vs 2 KO replicate?

I'm giving a sample of the output.

Any help or suggestion would be really appreciated.

```
dput(head(CTRL_PTPEM_Pol_II_pausing))
structure(list(`Transcript/RepeatID (cmd=analyzeRepeats.pl rna hg38 -strand both -count pausing -condenseGenes -d Control_Pol_II_new_set/ Control_Pol_II_set2/ PTPN6_KO_Pol_II_PTPEN/ PTPN6_KO_Pol_II_set2/ -rpkm)` = c("NM_173803",
"NM_001301039", "NR_132775", "NM_022114", "NM_001011666", "NR_160728"
), chr = c("chr16", "chr1", "chr3", "chr1", "chr7", "chr19"),
    start = c(15395704, 150272234, 13618331, 3069153, 28686083,
    22527485), end = c(15400754, 150277284, 13623381, 3074203,
    28691133, 22532535), strand = c("+", "+", "+", "+", "+",
    "-"), Length = c(250, 250, 250, 250, 250, 250), Copies = c(2,
    1, 1, 1, 1, 1), `Annotation/Divergence` = c("MPV17L|M-LPH|MLPH1|MLPH2|MPV17L1|-|16p13.11|protein-coding",
    "C1orf54|-|-|1q21.2|protein-coding", "SNORA93|-|-|3p25.1|snoRNA",
    "PRDM16|CMD1LL|KMT8F|LVNC8|MEL1|PFM13|-|1p36.32|protein-coding",
    "CREB5|CRE-BPA|CREB-5|CREBPA|-|7p15.1|protein-coding", "LOC105376917|-|-|19p12|ncRNA"
    ), `Control_Pol_II_new_set/ Pausing Ratio` = c(0.982, 1.832,
    1.114, 1.247, 3.476, 0.632), `Control_Pol_II_new_set/ Promoter Reads (-50,200)` = c(0.125,
    0.748, 0.249, 0.873, 0.499, 0), `Control_Pol_II_new_set/ GeneBody Reads (200,5000)` = c(0.13,
    0.338, 0.208, 0.669, 0.032, 0.091), `Control_Pol_II_set2/ Pausing Ratio` = c(0.686,
    1.131, 1.054, 0.816, 2.205, 1.005), `Control_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.1,
    0.7, 0.1, 0.3, 0.3, 0.1), `Control_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.203,
    0.604, 0.089, 0.396, 0.068, 0.099), `PTPN6_KO_Pol_II_PTPEN/ Pausing Ratio` = c(1.08,
    1.362, 1.108, 1.224, 1.943, 1.2), `PTPN6_KO_Pol_II_PTPEN/ Promoter Reads (-50,200)` = c(0.151,
    0.602, 0.151, 0.301, 0.452, 0.151), `PTPN6_KO_Pol_II_PTPEN/ GeneBody Reads (200,5000)` = c(0.125,
    0.392, 0.118, 0.212, 0.141, 0.094), `PTPN6_KO_Pol_II_set2/ Pausing Ratio` = c(1.652,
    0.909, 1.551, 1.395, 1.774, 0.584), `PTPN6_KO_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.429,
    0.536, 0.429, 0.536, 0.321, 0.107), `PTPN6_KO_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.206,
    0.603, 0.229, 0.346, 0.123, 0.279)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

```

# Answer
> 1 votes
but I would say the pausing index is something that is invented, commonly used, but with no clear indication of what is the expected distribution or how it can be subjected to test like you've mentioned.

I assume you would like to get your work published, so going by the papers with IF \> 5 such as this, and this, and this, the clear strategy is use a boxplot to show the difference between conditions, and apply a statistical test such as mann whitney u or t.test to say, there is a difference.

Let's try this, sorry not very good in R,using your dataframe which i named x:

```
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

df = x %>% select(contains("pausing") | contains("Transcript")) %>%
rename_at(vars(1),substr,1,10) %>%
pivot_longer(-Transcript,names_to = "condition") %>%
mutate(condition=case_when(
grepl("Control",condition) ~ "Control",
grepl("PTPN6",condition) ~ "PTPN6_KO"
)) %>%
group_by(Transcript,condition) %>% summarize_all(mean)

```

In the data frame above, I have averaged the pausing index between replicates (not very good but then again, most papers have only 1 replicate), and group them by condition.

```
# A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   <chr>        <chr>     <dbl>
 1 NM_001011666 Control   2.84
 2 NM_001011666 PTPN6_KO  1.86
 3 NM_001301039 Control   1.48
 4 NM_001301039 PTPN6_KO  1.14
 5 NM_022114    Control   1.03
 6 NM_022114    PTPN6_KO  1.31
 7 NM_173803    Control   0.834
 8 NM_173803    PTPN6_KO  1.37
 9 NR_132775    Control   1.08
10 NR_132775    PTPN6_KO  1.33
11 NR_160728    Control   0.818
12 NR_160728    PTPN6_KO  0.892

```

Then we plot:

```
p <- ggboxplot(df, x = "condition", y = "value",
                color = "condition", palette =c("#c72e29","#016392","#be9c2e","#098154"),
                add = "jitter")
p+stat_compare_means(comparisons=list( c("Control", "PTPN6_KO")))

```

---
Tags: homer
---


title: Coronavirus phylogeny and evolution

text:
# Question
Title: Coronavirus phylogeny and evolution
This new article presents the phylogeny of the coronaviruses and the placement of the new coronavirus in it. This made me wonder of what we know about the coronavirus evolution, particularly where it comes to evading human immune defenses: it is not mutating as fast as HIV - indeed, the genomes responsible for outbreaks in various countries are almost identical. In the same time it is unlikely to undergo gene reshuffling typical of influenza, since coronaviruses have single RNA genome.

# Answer
> 1 votes
FYI coronaviruses have proof reading activity, unlike HIV, its a huge genome, unlike HIV. Thus the point mutation rate is unlikely to be as rapid as HIV.

The phylogenies in this paper are not really about the evolution or epidemiology of the virus. The trees used for viral taxonomy and form part of the International Committee for the Taxonomy of Viruses (ICTV) assessment in officially naming the virus. The ICTV is the authority in naming a virus. The ICTV essentially called it SARS2 as opposed to the WHO who used the name COVID-19. Thats the point of this paper and will be massively cited.

The broader evolutionary epidemiology of "SARS2" is for further studies. The Lancet paper on 22nd February gave some indication of this (and I suspect you know about that paper).

---
Tags: phylogenetics, genome, phylogeny, covid-19, sars-cov-2
---


# Question
Title: Difference between clinical significance and statistical significance
I am working on a clinical data analysis project and am from IT domain.

While I am analyzing the data, I came to know that though some of the variables are not statistically significant, they are clinically significant.

I also read that difference in median values between case and control group for variables is used as a measure to indicate clinical significance

Can someone help me understand this with a simple example as to why difference in median between groups is considered for clinical significance?

# Answer
> 2 votes
Another viewpoint could be : statistically significant implies use of sampling theory to establish presence of significant difference between case and control groups. The "clinical " significance seems to test for correctness or validity of measures of variables. Apparenty, the measurement theory is meant for checking validity of data - the difference could arise on account of incorrect measurements. The non-parameteric tests - it seems - take care of measurement errors and hence, produce perhaps more reliable results.

# Answer
> 3 votes
This is non-parametric statistics, the mean requires conformation to the normal distribution (but there are exceptions), or perhaps better put it requires a fixed relationship between mean and standard deviation. Non-parametric stats always use the median and are distribution free. Wilcoxon test/ Mann-Whitney U test are used to replace t-test. The Kruskal-Wallis test replaces the 1-way ANOVA Finally the Friedman test replaces 2-way ANOVA but has its critiques. The t- and ANOVA tests require normality, F-max is one approach to determine normality.

---
Tags: database, errors, validation, p-values
---


# Question
Title: Prediction model for spread of coronavirus
A previous thread on Bioinformatics StackExchange citing an article by the Economist which references a cartoon of predictive outcomes by the CDC (Economist article on coronavirus) suggests a model in epidemiology using the "basic reproductive rate", commonly known as R0. Although I'm no bioinformatician, I imagine that one could tweak a model based on similar infections like SARS and HDN1. Yet, I found no models predicting or giving confidence intervals for the spread of coronavirus. I wonder if the the data is still too sparse, e.g. we don't know if it's seasonal (mentioned in the above thread) and we cannot infer the virulence.

Would such a model be reliable? If not, how many months of data would we need to reliably use a model predicting the spread?

# Answer
> 2 votes
The original R0 calculation was developed in the early 1900s to model malaria transmission and is known as the Ross-MacDonald model so your idea is over 100 years old. The model was very famous for identifying the optimal strategy for malarial control by targeting the adult mosquito rather than the mosquito larvae. What it said was reducing the mosquito population and reducing malaria transmission were two different issues. So the bednet strategy is based in this idea and particularly spraying the inside of houses with insecticide.

R0 was famously used in the 1980s to predict the HIV pandemic when WHO thought it would remain localised in Africa.

The problem with R0 and Coronaviruses is we don't know the asymptomatic population who could also transmit the virus.

This population is potentially very large. Epidemiologists will attempt extrapolation from cruise ship quarrentine, eg the Diamond Princess where everyone is tested to understand this transmission parameter and project this onto the rest of the planet. Clearly this could have limitations or inaccuracies. The other issue I am baffled by is seasonal transmission and whether this will occur with COVID-19.

One thing R0 of COVID-19 will identify is the latent period transmission which occurs before the symptoms, i.e. the incubation period. Essentially, people with cough/mild symptoms before they self isolate and become really ill with a lung infection will transmit the virus. In this context R0 will be/was very influential for understanding when to move between quarrentine and containment phases of the pandemic.

---
Tags: statistics, covid-19, sars-cov-2, coronavirus, epidemiology
---


# Question
Title: Tools or softwares for clustering evaluation
I am using a few clustering algorithms as well as my own tool for protein sequences. Also, I have benchmarked clustering results (i.e. ground truth) to compare with the results. Is there any tool/software/script to generate evaluation metrics for the clustering results.

The clustering results/benchmarked are in the following format

```
1  Name1
1  Name2
2  Name3
2  Name4
2  Name5
.....

```

UPDATE: I am evaluating the clustering performance using scikit

# Answer
The definitive answer is no, that is why they are not so popular. There are metrics that can be used:

* you can examine "minimum distance" as a criterion under the assumption that the lowest distance between taxa represents the largest homology (e.g. if it is sequence data), which would be translate to other biological scenarios.
* resampling, e.g. bootstrapping, jack-knifing and examine the tree structure of each node against the biological plausibility of the answer being correct.

For both points a biological a priori, i.e. biological plausibility of assessing the tree structure is required to confirm/reject the clustering method.

Clustering algorithms are not like the probabilistic methods, where the highest likelihood is the singular criteria and of course they do give different answers. Its a brief explanation I know but hope it makes sense.
> 1 votes

# Answer
I am evaluating the clustering performance using scikit
> 1 votes

---
Tags: sequence-analysis, clustering
---


# Question
Title: How do I Investigate and test the hypothesis on the effect of street-light regime on insects?
I need advice on a suitable analysis of a dataset that has an akward design because we investigate a real-world situation. We have two sites (A, B) where we can conduct the research. The hypothesis is: under-dimmed light, we catch less number of insects than under-fully lit street lights.

The two sites each have ten street lights that were sampled during 16 nights of dimmed light and during 16 nights of full light. During these total of 32 nights, we sampled insects at each street-light pole (10 per site, 20 in total). For A and also at B, we also have temperature as a variable (nightly average, max, min) - a measure per site not per light pole.

I would like to fit a model: number of insects ~ dimming, temp\_mean, temp\_min, temp\_max, site

I include "site" as categorical variable and would like to use the 32 nights (16 dimmed, 16 fully lit) as 32 unique samples (each night as one sample). I am not sure whether this is "statistically sound". If not, does anybody have an idea on how to evaluate such data in order to prove the hypothesis? Would you go for t-test, PCA or something else?

```
streetlight <- data.frame(site=c(rep("A",16),rep("B",16),rep("A",16),rep("B",16)),
                          dimming=c(rep(1,32),rep(0,32)),temp_mean=runif(32, min=20, max=24),t
                          emp_min=runif(32, min=12, max=19),temp_max=runif(32, min=25, max=28),
                          insects=round(runif(32, min=5, max=20),0))

```

# Answer
> 1 votes
The hypothesis can be examined by t-test. You can group the data (number of catches) into dimmed light, and fully lit street lights and proceed with t-test to assess whether dimmed light results in less en-catchment or higher en-catchment or equivalent encatchment. Alternatively, a linear model can be applied subject to certain conditions here. PCA is not appropriate here.

# Answer
> 0 votes
Please avoid PCA. Go with logistic regression and work out the odds ratio. This is far and away the easiest most powerful method.

Oh and you're likely to see a negative binomial distribution, very common for insect distributions. Your next step would be to get some DNA out of them

---
Tags: r, linear-regression
---


# Question
Title: Identifying proteins in a concatenated protein sequences
I have a 120-concatenated protein sequence in a fasta file. I would like to split by proteins. How can I identify each protein?

# Answer
> 1 votes
The only way I can think is to break the protein apart on methionine residues, blast each amino acid locus and concatenate any locus that does not score a query protein position of '1'. The theory is simple every protein starts with a methionine residue, however there are plenty of M residues within the protein that don't represent the beginning of the gene.

Breaking the protein apart on M would be easy, but then automating and parsing the blast would be fiddly. Obvioulsy you will need to code it and if its a one off would it be worth it?

---
Tags: sequence-alignment, fasta, proteins
---


# Question
Title: What is a proper way for random subsampling of metagenomic data?
Let's say we have a metagenomic sample that is paired-end FASTQ files including 10,000,000 DNA reads collected using shotgun sequencing. How would one make a random subsample of the mentioned metagenomic sample with for example 1,000,000 reads? I know about the `seqtk` function that could be used in such cases, but my question is more about whether it is meaningful to represent a metagenomic sample with a small subsample generated using this way or it has to satisfy certain criteria, e.g. in terms of the number of reads and diversity to be a valid representative of the metagenomic sample? Knowing that some species might be very rare in the sample, would it ever be proper to do subsampling on metagenomic samples?

# Answer
I am not NGS expert, I known ML. Essentially you are performing data augmentation which is essential for ML. What I question is whether a replicate of e6 from a population of e7 would really be sufficient for ML because you need sample sizes (replicates in this case) of \>100000.

The answer is really simple you bootstrap the data (sample with replacement), to obtain your replicates. Numpy has a random module that will do this,

```
p.random.seed(345)
pop = np.random.randint(0,10000000 , size=1000000)

```

Use `pop` to select from fastaq. The other way in is to define each fastq sequence within pop then simply,

```
sample = np.random.choice(pop, size=1000000)

```

Or a combination of the two. If you can do ML you can do Python coding.
> 1 votes

---
Tags: fastq, taxonomy, microbiology
---


# Question
Title: Prediction of prokaryotic origins of replication (ORI)
I want to predict origins of replication (ORI) on hundreds of prokaryotic genomes. The most straight-forward solution would be to use most commonly used tool, Ori-Finder.

It uses integrated gene prediction, analysis of base composition asymmetry, distribution of DnaA boxes, occurrence of genes frequently close to oriC regions and phylogenetic relationships to predict ORIs.

However, I came across a paper by Parikh et al which describes a single measure (correlated entropy measure / CEM) with which they achieved 100% ORI-identification on 500 genomes.

I like the CEM approach more than a close-source Windows-only black-box.

Parikh tested their algorithm on NCBI genomes. How are ORI annotated on these reference genomes? I couldn't find such an annotation. Does the sequence simply begin at the ORI in leading strand direction by convention? Were the ORI of the reference genomes identified by experiment?

# Answer
> 1 votes
I have not seen the papers but I know the technique SVM. This is a machine learning approach its pretty good, but has been around this the 1980s.

Its quite simple known ORIs are given with the associated genetic trait, e.g. GC skew, SVM will then regress the GC skew of a known ORI against the GC content of the genome, there will be a training, parameter optimization and test split in the data, typically 80:20:20. De novo data will then be feed into the trained algorithm, in the same way the test data was, and a prediction will be generated. There are weaknesses, the authors have used only one ML technique, typically numerous will be used e.g. SVM, random forest, lasso regression etc... In addition, the prediction accuracy for GC is low (75%), in addition 500 genomes isn't great for ML, you need 1000s.

---
Tags: genome, algorithms
---


# Question
Title: Removing common variants in the 1000 genomes database from .vcf
I have 15 .vcf files. I need to remove \`common variants in the 1000 genomes database' appearing in at least 0.5% of the population

Do you know from where I may start?

Thank you so much

# Answer
Use VEP, ANNOVAR or snpEff to annotate your VCF file (I'd recommend combining your VCFs into a single file if they're single sample VCFs or are all comprised of samples from the same experiment/cohort). You can annotate the VCF file with 1000g (among a ton of other annotation sources). Once done, you can use `bcftools view` to subset the VCF as required.

Or, you can subset the standalone 1000g VCF and pick regions with your frequency threshold, then you can use that regions file to subset your VCF file(s).

Remember to left-align and normalize (and split multi-allelic records in) your VCF files so all comparisons can be done accounting for chr, pos, ref and alt fields.
> 4 votes

# Answer
Think removing everything \< 0.05 from `1000g2015aug_all` column of ANNOVAR annotation file helps with this
> -3 votes

---
Tags: vcf, wgs, 1000genomes
---


title: Installing a package on Mac

text:
# Question
Title: Installing a package on Mac
I am trying to instll sigproextractor package on my Mac

I tried

```
pc-133-235:~ fi1d18$ python --version
Python 2.7.10
pc-133-235:~ fi1d18$ python3 --version
Python 3.8.2
pc-133-235:~ fi1d18$ pip --version
-bash: pip: command not found
pc-133-235:~ fi1d18$ curl https://bootstrap.pypa.io/get-pip.py -o get-pip.py

```

But whatever I am trying says that can't find pip

Any help?

# Answer
> 4 votes
maybe you could try `python3 -m pip install package`

# Answer
> 2 votes
At a guess you don't have python3 installed and hence pip3

https://stackoverflow.com/questions/17271319/how-do-i-install-pip-on-macos-or-os-x

Personally I would install via Anaconda3 and skip a direct Python3 installation altogether. The answer above recommends using homebrew, again personally I wouldn't touch it and go straight for a full Anaconda3 installation. You might be able to install via

```
conda install -c XXX sigproextractor

```

Where XXX is the channel, but its certainly on github. Anaconda3 is very easy to install.

---
Tags: python, conda
---


# Question
Title: How to specify resources for cluster in snakemake
I am trying to run `snakemake` version `3.5.4` on my cluster, where jobs are partially defined by a common template build according manual and also defined though rule-specific parameters. But for some reason cluster configuration is not expanding wildcards (and command does). When I run

```
snakemake data/genome.fa.gz.bwt -p --jobs 10 --cluster-config cluster.json --cluster "bsub -J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} -o {cluster.output} -e {cluster.error}"

```

`-J "{rule}.{wildcards}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host}` expands as expected. But `-o {cluster.output} -e {cluster.error}` is replaced by literally `logs/cluster/{rule}.{wildcards}.out` and `logs/cluster/{rule}.{wildcards}.err`, i.e. there is no expansion in config file (the config file is attached bellow).

I am puzzled since exactly this lines are in the documentation. What am I doing wrong? Why `snakemake` is not expanding wildcards in the config file?

Here is my `cluster.json`

```
{
    "__default__" :
    {
        "queue"    : "bgee",
        "one_host" : "\"span[hosts=1]\"",
        "name"      : "{rule}.{wildcards}",
        "output"    : "logs/cluster/{rule}.{wildcards}.out",
        "error"     : "logs/cluster/{rule}.{wildcards}.err"
    },
}

```

Here is relevant part of `Snakefile` :

```
rule index_reference :
        threads: 1
        resources: M=20000000
        input : "data/{reference}/genome.fa.gz"
        output : "data/{reference}/genome.fa.gz.bwt"
        shell : "scripts/index_genome.sh {input}"

```

# Answer
> 5 votes
This feature was added in version 3.6.0. Look at release notes.

`{wildcards}` normally expands to a something that looks vaguely like a dictionary. In other words, something like `key1=value1,key2=value2`. Whether you really want to use that as a log file name is up to you.

Let's use a very contrived example to demonstrate this. Here is the `Snakefile`:

```
rule all:
    input: "input1.human.txt"

rule blah:
    output: "{sample}.{reference}.txt"
    shell: "touch {output}"

```

I'm using your `cluster.json`, so with that let's use a contrived cluster submission program:

```
snakemake --cluster-config cluster.json --cluster "echo 'rule.wildcards: {rule}.{wildcards} cluster: {cluster} -o {cluster.output} -e {cluster.error} -n {cluster.name}' >> log && " --jobs 1

```

If you run that and look in `log`, you can see that `{rule}` is expanded as expected, `{wildcards}` is expanded to `sample=input1,reference=human`. The `output` and `error` parts from `cluster.json` are expanded correctly for me. Below are the contents of `log`:

```
rule.wildcards: blah.sample=input1,reference=human cluster: one_host="span[hosts=1]",output=logs/cluster/blah.sample=input1,reference=human.out,name=blah.sample=input1,reference=human,queue=bgee,error=logs/cluster/blah.sample=input1,reference=human.err -o logs/cluster/blah.sample=input1,reference=human.out -e logs/cluster/blah.sample=input1,reference=human.err -n blah.sample=input1,reference=human

```

I'm not a big fan of using `{wildcards}` in paths unless you're doing a parameter sweep, but to each their own. Note that the above results are from snakemake-4.3.1.

# Answer
> 3 votes
The solutions are either, upgrade `snakemake` or just take all wildcards out of `cluster.json` file and write them directly to snakemake execution command :

```
snakemake -p --jobs 10 --cluster-config cluster.json --cluster "bsub -J "{rule}" -q {cluster.queue} -n {threads} -M {resources.M} -R {cluster.one_host} -o {"logs/cluster/{rule}.{wildcards}.out" -e "logs/cluster/{rule}.{wildcards}.err""

```

However, I ended up instead of using of json config file in writing a simple bash wrapper `snakemake_cluster.sh` that remembers the execution command for me in more readable format :

```
#!/usr/bin/env bash

snakemake $@ -p --jobs 10 --cluster "bsub \
-J {rule} \
-q bgee \
-n {threads} \
-M {resources.M} \
-R \"span[hosts=1]\" \
-o logs/cluster/{rule}.{wildcards}.out \
-e logs/cluster/{rule}.{wildcards}.err"

```

Edit : All spaces after `--cluster` command got to be escaped, otherwise it would add arguments to `snakemake` not to `bsub`.

Edit 2 : `$@` means all arguments, therefore it can be either targets, flags or combinations (targets got to be specified first). So, `snakemake_cluster.sh my_file` builds `my_file`, `snakemake_cluster.sh --quiet` builds default target quietly, and `snakemake_cluster.sh my_file --quiet` builds `my_file` quietly.

# Answer
> 1 votes
As of the most recent snakemake (v5.10.0) the recommended way of defining resources is through a profile.

It seems you are trying to submit a job on an LSF cluster system - there is a snakemake profile for LSF to save you having to do this work yourself.

It reduces your snakemake submission to something like

```
snakemake --profile "lsf"

```

---
Tags: snakemake
---


# Question
Title: polishing assembled genome to QV50 value
I am trying to assemble some genomes from pacbio reads using HGAP3 followed by quiver polishing. However, even after multiple rounds of quiver polishing, I cannot attain quality value 50 (QV50) for a genome.

QV50 seems to be recommended and I am able to get upto QV30 only. I have done polishing for almost 15 times already. I have tried polishing with different parameter changes like subread length cutoff and blasr maximum divergence.

Can anyone suggest me on how I can achieve an assembled genome of QV50 value?

Thank you

# Answer
> 1 votes
15 times is too many times. If you were going to get quiver to get you good QVs it would have happened after 2-3 rounds probably.

Some questions:

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

---
Tags: assembly, pacbio
---


# Question
Title: Translating a genome sequence into possible frames
I have a sequence below from MYC gene about which I need to translate the sequence in all possible frames AND identify (for each frame) which codons are actually used in MYC

```
AAAGCAATAATACAATTTAAAACCTGGGTCTCTAGAGGTGTTAGGACGTG
GTGTTGGGTAGGCGCAGGCAGGGGAAAAGGGAGGCGAGGATGTGTCCGAT
TCTCCTGGAATCGTTGACTTGGAAAAACCAGGGCGAATCTCCGCACCCAG
CCCTGACTCCCCTGCCGCGGCCGCCCTCGGGTGTCCTCGCGCCCGAGATG
CGGAGGAACTGCGAGGAGCGGGGCTCTGGGCGGTTCCAGAACAGCTGCTA
CCCTTGGTGGGGTGGCTCCGGGGGAGGTATCGCAGCGGGGTCTCTGGCGC
AGTTGCATCTCCGTATTGAGTGCGAAGGGAGGTGCCCCTATTATTATTTG
ACACCCCCCTTGTATTTATGGAGGGGTGTTAAAGCCCGCGGCTGAGCTCG
CCACTCCAGCCGGCGAGAGAAAGAAGAAAAGCTGGCAAAAGGAGTGTTGG
ACGGGGGCGGTACTGGGGGTGGGGACGGGGGCGGTGGAGAGGGAAGGTTG
GGAGGGGCTGCGGTGCCGGCGGGGGTAGGAGAGCGGCTAGGGCGCGAGTG
GGAACAGCCGCAGCGGAGGGGCCCCGGCGCGGAGCGGGGTTCACGCAGCC
GCTAGCGCCCAGGCGCCTCTCGCCTTCTCCTTCAGGTGGCGCAAAACTTT
GTGCCTTGGATTTTGGCAAATTGTTTTCCTCACCGCCACCTCCCGCGGCT
TCTTAAGGGCGCCAGGGCCGATTTCGATTCCTCTGCCGCTGCGGGGCCGA
CTCCCGGGCTTTGCGCTCCGGGCTCCCGGGGGAGCGGGGGCTCGGCGGGC
ACCAAGCCGCTGGTTCACTAAGTGCGTCTCCGAGATAGCAGGGGACTGTC
CAAAGGGGGTGAAAGGGTGCTCCCTTTATTCCCCCACCAAGACCACCCAG
CCGCTTTAGGGGATAGCTCTGCAAGGGGAGAGGTTCGGGACTGTGGCGCG
CACTGCGCGCTGCGCCAGGTTTCCGCACCAAGACCCCTTTAACTCAAGAC
TGCCTCCCGCTTTGTGTGCCCCGCTCCAGCAGCCTCCCGCGACGATGCCC
CTCAACGTTAGCTTCACCAACAGGAACTATGACCTCGACTACGACTCGGT
GCAGCCGTATTTCTACTGCGACGAGGAGGAGAACTTCTACCAGCAGCAGC
AGCAGAGCGAGCTGCAGCCCCCGGCGCCCAGCGAGGATATCTGGAAGAAA
TTCGAGCTGCTGCCCACCCCGCCCCTGTCCCCTAGCCGCCGCTCCGGGCT
CTGCTCGCCCTCCTACGTTGCGGTCACACCCTTCTCCCTTCGGGGAGACA

```

Where should I look for all 4 pieces of information easily? Please guide.

# Answer
> 1 votes
The 6-frame translation results are below you will need to do the rest. It is translated 5' to 3' but appears to contain two splice sites. You will need to use blast to resolve the start/stop of the splicing.

All 6-codon translation is available the Expasy server here. You will need to blast to resolve the splicing.

---

5'3' Frame 1 KAIIQFKTWVSRGVRTWCWVGAGRGKGRRGCVRFSWNR-LGKTRANLRTQP-LPCRGRPRVSSRPRCGGTARSGALGGSRTAATLGGVAPGEVSQRGLWRSCISVLSAKGGAPIII-HPPCIYGGVLKPAAELATPAGERKKKSWQKECWTGAVLGVGTGAVEREGWEGLRCRRG-ESG-GASGNSRSGGAPARSGVHAAASAQAPLAFSFRWRKTLCLGFWQIVFLTATSRGFLRAPGPISIPLPLRGRLPGFALRAPGGAGARRAPSRWFTKCVSEIAGDCPKGVKGCSLYSPTKTTQPL-GIALQGERFGTVARTARCARFPHQDPFNSRLPPALCAPLQQPPATMPLNVSFTNRNYDLDYDSVQPYFYCDEEENFYQQQQQSELQPPAPSEDIWKKFELLPTPPLSPSRRSGLCSPSYVAVTPFSLRGD 5'3' Frame 2 KQ-YNLKPGSLEVLGRGVG-AQAGEKGGEDVSDSPGIVDLEKPGRISAPSPDSPAAAALGCPRARDAEELRGAGLWAVPEQLLPLVGWLRGRYRSGVSGAVASPY-VRREVPLLLFDTPLVFMEGC-SPRLSSPLQPARERRKAGKRSVGRGRYWGWGRGRWRGKVGRGCGAGGGRRAARARVGTAAAEGPRRGAGFTQPLAPRRLSPSPSGGAKLCALDFGKLFSSPPPPAAS-GRQGRFRFLCRCGADSRALRSGLPGERGLGGHQAAGSLSASPR-QGTVQRG-KGAPFIPPPRPPSRFRG-LCKGRGSGLWRALRAAPGFRTKTPLTQDCLPLCVPRSSSLPRRCPSTLASPTGTMTSTTTRCSRISTATRRRTSTSSSSRASCSPRRPARISGRNSSCCPPRPCPLAAAPGSARPPTLRSHPSPFGET 5'3' Frame 3 SNNTI-NLGL-RC-DVVLGRRRQGKREARMCPILLESLTWKNQGESPHPALTPLPRPPSGVLAPEMRRNCEERGSGRFQNSCYPWWGGSGGGIAAGSLAQLHLRIECEGRCPYYYLTPPLYLWRGVKARG-ARHSSRREKEEKLAKGVLDGGGTGGGDGGGGEGRLGGAAVPAGVGERLGREWEQPQRRGPGAERGSRSR-RPGASRLLLQVAQNFVPWILANCFPHRHLPRLLKGARADFDSSAAAGPTPGLCAPGSRGSGGSAGTKPLVH-VRLRDSRGLSKGGERVLPLFPHQDHPAALGDSSARGEVRDCGAHCALRQVSAPRPL-LKTASRFVCPAPAASRDDAPQR-LHQQEL-PRLRLGAAVFLLRRGGELLPAAAAERAAAPGAQRGYLEEIRAAAHPAPVP-PPLRALLALLRCGHTLLPSGR 3'5' Frame 1 CLPEGRRV-PQRRRASRARSGG-GTGAGWAAARISSRYPRWAPGAAARSAAAAGRSSPPRRSRNTAAPSRSRGHSSCW-S-R-GASSREAAGAGHTKREAVLS-RGLGAETWRSAQCAPQSRTSPLAELSPKAAGWSWWGNKGSTLSPPLDSPLLSRRRT--TSGLVPAEPPLPREPGAQSPGVGPAAAEESKSALAPLRSRGRWR-GKQFAKIQGTKFCAT-RRRREAPGR-RLREPRSAPGPLRCGCSHSRPSRSPTPAGTAAPPNLPSPPPPSPPPVPPPSNTPFASFSSFSRRLEWRAQPRALTPLHKYKGGVK---GHLPSHSIRRCNCARDPAAIPPPEPPHQG-QLFWNRPEPRSSQFLRISGARTPEGGRGRGVRAGCGDSPWFFQVNDSRRIGHILASLFPCLRLPNTTS-HL-RPRF-IVLLL 3'5' Frame 2 VSPKGEGCDRNVGGRAEPGAAARGQGRGGQQLEFLPDILAGRRGLQLALLLLLVEVLLLVAVEIRLHRVVVEVIVPVGEANVEGHRRGRLLERGTQSGRQS-VKGVLVRKPGAARSARHSPEPLPLQSYPLKRLGGLGGGIKGAPFHPLWTVPCYLGDALSEPAAWCPPSPRSPGSPERKARESAPQRQRNRNRPWRP-EAAGGGGEENNLPKSKAQSFAPPEGEGERRLGASGCVNPAPRRGPSAAAVPTRALAALLPPPAPQPLPTFPLHRPRPHPQYRPRPTLLLPAFLLSLAGWSGELSRGL-HPSINTRGVSNNNRGTSLRTQYGDATAPETPLRYLPRSHPTKGSSCSGTAQSPAPRSSSASRARGHPRAAAAGESGLGAEIRPGFSKSTIPGESDTSSPPFSPACAYPTPRPNTSRDPGFKLYYCF 3'5' Frame 3 SPRREKGVTAT-EGEQSPERRLGDRGGVGSSSNFFQISSLGAGGCSSLCCCCW-KFSSSSQ-KYGCTES-SRS-FLLVKLTLRGIVAGGCWSGAHKAGGSLELKGSWCGNLAQRAVRATVPNLSPCRAIP-SGWVVLVGE-REHPFTPFGQSPAISETHLVNQRLGARRAPAPPGARSAKPGSRPRSGRGIEIGPGALKKPREVAVRKTICQNPRHKVLRHLKEKARGAWALAAA-TPLRAGAPPLRLFPLAP-PLSYPRRHRSPSQPSLSTAPVPTPSTAPVQHSFCQLFFFLSPAGVASSAAGFNTPP-IQGGCQIIIGAPPFALNTEMQLRQRPRCDTSPGATPPRVAAVLEPPRAPLLAVPPHLGREDTRGRPRQGSQGWVRRFALVFPSQRFQENRTHPRLPFPLPAPTQHHVLTPLETQVLNCIIA

---
Tags: sequence-alignment, ngs, proteins, blast, phylogenetics
---


# Question
Title: How to construct the phylogenetic tree of the "High Recombination" genus streptomyces?
I need to construct the phylogenetic tree of the whole genus of streptomyces. And **my goal** is to **find out how my target strep species distributed in the whole strep genus and obtain the distribution pattern**. Strep is a big genus, has 658 species(65% of them have the genome info), and my target species is about 25. **But the traditional 16S method seems not helping when it comes to strep is a genetic recombination genus**. And still lack of MLST scheme to distinguish my target species. I want to get the phylogeny for the whole genus, but I lack related genetics knowledge.

There is something I did before try to get the tree.

1. I ran the orthomcl program for 400 strep genomes which I got from NCBI, and the result is: I got 22 orthologues
2. I ran the fastani program for 200 strep genomes(type material) and got the ani result.

What can I do to achieve my goal? Thanks!

# Answer
> 1 votes
You need bootstrap support using a model-based tree building algorhithm, via maximum likelihood (a few people use Bayes). The file format is relaxed phylip format (please submit a separate question if you have difficulties here - its a bit tricky).

I use standard RAxML here, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

```
raxmlHPC -m GTRGAMMA -n output.tre -o outgroup-p 10000 -s alignment_file.phy

```

To make a bootstrap tree,

```
raxmlHPC -b 459 -m GTRGAMMA -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

```

This tree will bootstrap for 500 replications, however to start I would use 100 replications.

Make a consensus tree of the bootstraps,

```
raxmlHPC -f b -m GTRGAMMA -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

```

You require bootstrap support \>80% and please repeat this with and without 5-2/5-3 (it still looks long)

The cluster you have access to is fine for the calculation, it will take aound 24 hours to complete one bootstrap calculation for one data set and obviously you need to parallelise you calculation across the 22 contigs.

Viewing the tree, FigTree (for Mac OSX) is easy.

**Root**

Rooting can be complicated because I don't really know your bacteria.

---

The recombination issue is more complicated, but I would construct 22 trees from your contigs and assess them for congruence. Panmixia is a concerned, which means too much recombination

---
Tags: fasta, phylogenetics, phylogeny, bacteria
---


# Question
Title: Cutting the sequence into several sequences with the information of a dataframe
I have a fasta file with several 120-concatenation protein sequences. I also have a data frame with 120 names of proteins in column 1 and their length in column 2. Using this dataframe I want to create a new fasta file per sequence in which I get each protein independent with the name of the protein and the bacteria.

Here my fasta file.

```
Bacteria_1
Seq('GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...IEE', SingleLetterAlphabet())

Bacteria_2
Seq('GRRKDSVARVYLRPGKGEIVVNKRPVEDYFGRETLKMVLRQPLELTESLDKFDI...YK-', SingleLetterAlphabet())

```

Here my dataframe

```
      Names  LEN
0     Ribosomal_S9:  121
1     Ribosomal_S8:  129
2    Ribosomal_L10:  100
3             GrpE:  166
4           DUF150:  141
..              ...  ...
115      TIGR03632:  117
116      TIGR03654:  175
117      TIGR03723:  314
118      TIGR03725:  212
119      TIGR03953:  188

[120 rows x 2 columns]

```

The name of the genes and their length is the same for all my bacteria so this df is applicable for all sequences in my fasta file.

I am trying to get new fasta files like this:

```
>Bacteria_1 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_1 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_1 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

```

The second fasta file

```
>Bacteria_2 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_2 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_2 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

```

# Answer
We wouldn't write the code for you from scratch (thats your job), but its an easy problem because in Python any string is an "array" or "dictionary". Thus

```
bacteria1 = 'GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDVIEEXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
print (bacteria1[0:120])

```

This will define the first gene in your list, remembering the array starts at 0 whilst your gene starts at position 1 (so you need to subtract 1). You simply parse the input list so the gene lengths form a dictionary and loop through the dictionary of gene lengths ensuring the last length is added to the current one and then issuing a write command.

You don't really need biopython in this instance, unless it will parse the gene length file you need to open (I'm not aware it can).

---

BTW the format you supplied isn't fasta it looks like a Biopython statement, again in this problem I wouldn't bother because Python is perfectly capable of splicing a string.
> 1 votes

---
Tags: fasta, ngs, python
---


# Question
Title: are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?
I have recently ran in R four geneset enrichment analysis in the same database (TCGA:breast cancer) comparing two intrinsic subtypes. The methods I used were:

* MIGSA, that imports mGSZ package and combines it with a SEA algorithm. Using RNA-seq TMM normalized counts as input, because the software supports it. All the parameters are as default.
* mGSZ with standard parameters, using RNA-seq TMM normalized log2 cpm. All parameters as default
* FGSEA preranked, using log fold change as ranking metric, obtained from limma deg from RNA-seq TMM normalized counts. All the parameters are as default.
* GSVA + limma: extracted by-sample enrichment score from RNA-seq TMM normalized log2 cpm, then ran a differential expression analysis with geneset enrichment scores as if they were simple gene expression values. I know this is not a formal approach but I wanted to compare. All the parameters are as default.

The genesets tested are the full combined list of KEGG, GO:BP, GO:MF and GO:CC. In all cases I trimmed the genesets to genesets of length \<500, which leaves a total of 21723 genesets tested.
Sadly, results are confusing. For starters mGSZ for some reason is cropping the output and giving results for only around 9700 of those genesets -I've already sent an email to the package mainteiner asking about that. But that isn't the biggest problem, the most shocking thing is that none of the analyses performedt has a good correlation to the results of any of the other analyses. Here I plotted the -log10(adjusted p-value + 0.0001) of each analysis agaisnt every other analysis, with lines indicating -log10(0.05)=1.3 and -log10(0.01)=2.
I was expecting the results to be somehow different, but here I can't see even something remotely close to a correlation between the different results. It is also evident that MIGSA and GSVA+limma output a lot more significatively enriched genesets than the other two methods.
I have gone through the code time and again, and I can't find anything errors. Also I have searched for studies doing this kind of direct comparison between methods, without luck so far. Has anybody experienced this kind of inconsistencies between gene set enrichment analysis results with different tools? What could I possible be doing wrong?.

# Answer
> 2 votes
Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.

Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.

This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)

The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): Toward a gold standard for benchmarking gene set enrichment analysis.

Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.

---
Tags: r, gsea, go-enrichment
---


title: Error in ruuning a package in R

text:
# Question
Title: Error in ruuning a package in R
I want to install `Signature Tools Lib R package`

I typed

```
> devtools::install()
Error: Could not find package root, is your working directory inside a package?
>

```

Packages does need this list of sotware

```
R package dependencies:
    VariantAnnotation,
    BSgenome.Hsapiens.UCSC.hg38,
    BSgenome.Hsapiens.1000genomes.hs37d5,
    BSgenome.Mmusculus.UCSC.mm10,
    BSgenome.Cfamiliaris.UCSC.canFam3,
    SummarizedExperiment,
    BiocGenerics,
    GenomeInfoDb,
    NMF,
    foreach,
    doParallel,
    doMC,
    lpSolve,
    ggplot2,
    methods,
    cluster,
    stats,
    NNLM,
    nnls,
    GenSA,
    gmp,
    plyr,
    RCircos,
    scales,
    GenomicRanges,
    IRanges,
    BSgenome

```

I have installed all of these one by one but again when trying things like the below, I get error

```
Error in tabToSNVcatalogue(subs = tmpSNVtab, genome.v = "hg19") :
  could not find function "tabToSNVcatalogue"
> devtools::test()
Error: Could not find package root, is your working directory inside a package?
>
> library(signature.tools.lib)
Error in library(signature.tools.lib) :
  there is no package called ‘signature.tools.lib’
>

```

# Answer
> 4 votes
If you do `?devtools::install`, you will see that the `install()` function from the `devtools` package, not surprisingly, requires a package to be specified:

```
install(pkg = ".", reload = TRUE, quick = FALSE, build = !quick,
  args = getOption("devtools.install.args"), quiet = FALSE,
  dependencies = NA, upgrade = "ask", build_vignettes = FALSE,
  keep_source = getOption("keep.source.pkgs"), force = FALSE, ...)

```

So it is absolutely normal that `devtools::install()` does not install anything. Just specify the packages(s) with the `pkg` argument.

And this is how you install `Signature Tools`:

```
install.packages("remotes")
remotes::install_github("Nik-Zainal-Group/signature.tools.lib")

```

---
Tags: r
---


# Question
Title: Comparing phylogenetic models with different datasets
I'm a linguist interested in phylogenetic tree inference using language data. I'm posting here because I'm using Bayesian phylogenetic methods in my work (probably using BEAST and/or RevBayes). For the purposes of this question, please accept the assumption that phylogenetic trees are generally applicable to languages (I'm aware this premise may be subject to dispute). My question is intended to be methodological and not necessarily domain specific. You could imagine I'm constructing biological species trees using morphological data.

I have two datasets: A set of binary traits ($D\_{b}$) and a set of continuous traits ($D\_{c}$).

I'd like to compare the results of phylogenetic tree inference with **(a)** binary data alone vs. **(b)** binary data and continuous data together. Model parameters are the same, except of course **(b)** requires some extra parameters to model the evolution of continuous traits. In other words, both the data and model of **(a)** is a subset of **(b)**.

I know this probably sounds like a strange research question. The reason it's of interest to me is that the binary data represent a data type that is well-established in historical linguistics and often used to infer language trees. The continuous data come from different language structures and haven't been used to infer language trees before. I want to know if adding this new kind of data can help aid tree inference significantly or if it has a neutral effect (or even just introduces more noise).

I'm familiar with phylogenetic model comparison using Bayes factors. But in the standard use case I'm familiar with, the dataset stays exactly the same and Bayes factors are used to evaluate the fit of different evolutionary models.

**My question: Is there a standard way of evaluating the effect of adding more data to a phylogenetic tree model? Can marginal likelihoods be compared still even though the data differs in each model?**

If not, would it make any sense to reframe the question as one of model comparison? I'm imagining, for example, running both **(a)** and **(b)** with identical datasets ($D\_b + D\_c$ for both). In **(a)**, there would be a trivial model of trait evolution that doesn't aid tree inference at all—something like a model where traits can jump from any value to any other value with equal probability(?). In **(b)**, there would be a more standard model of continuous trait evolution, e.g. Brownian motion.

# Answer
The key issue is the data is subject to reticulation, languages are seldom clean breaks. In nucleotide data this can be identified due to deviations from physical linkage and call it recombination. Linkage doesn't apply to languages.

The fundamental assumption is one tree - the same tree - describes all the data, with reticulation that may not be a good assumption. Thus English is a hybrid language, and lots of English words find their way back into French and German 'parents'. Or perhaps a better example is Gaelic mingles with English over time. There are methods used to monitor reticulation in languages.

Thus id keep it simple. The branch length is not only directly comparable it provides linguistical information that is useful to understand. This is the amount of 'evolutionaey' change per word (it that is what you binary data is describing) and the relationship to the node ('breakpoiint') is of primary interest.

The second step is to simply draw a tree of data set A and compare it to data set A+B and data set B and the probabilities describing the sister groups. Are the trees congruent, i.e. are they statistically the same?

The area you describing is more given to parametric bootstrapping, using custom resampling, but you need to be certain that the model being used is appropriate for the data.

The final point is whether a reversible or directional matrix is being used, which impacts on the specific result.
> 1 votes

---
Tags: phylogenetics, statistics, phylogeny, modelling, beast
---


# Question
Title: What is the per-gene score in BUSCO output?
The BUSCO `_full_table.tsv` output has scores per gene (an example bellow).

In the original publication, authors just mention a condition to satisfy a bitscore, but I could not find what the bitscore actually is. In the use guide they say

> If found to be complete, whether single-copy or duplicated, the BUSCO matches have scored within the expected range of scores ...

but again, I don't understand what the score actually represents. I am even not sure if greater the better or visa reverse. Does anybody know what the score is?

```
# BUSCO version is: 4.0.2
# The lineage dataset is: insecta_odb10 (Creation date: 2019-11-20, number of species: 75, number of BUSCOs: 1367)
# Busco id      Status  Sequence        Gene Start      Gene End        Score   Length
47at50557       Complete        scf1:20448-42284    20448   42284   4607.1  3554    https://www.orthodb.org/v10?query=47at50557     Immunoglobulin-like domain
547at50557      Complete        scf2:9780-23590    9780    23590   2743.8  1555    https://www.orthodb.org/v10?query=547at50557    protein unc-80 homolog
755at50557      Complete        scf3:9036-28083    9036    28083   3149.9  2206    https://www.orthodb.org/v10?query=755at50557    SRCR-like domain

```

# Answer
You are correct, it is the score of the candidate protein against the BUSCO profile hmm assessed by hmmer. If it passes a threshold, it is validated as being that BUSCO gene. The threshold is defined for each profile based on the score of the species used to produce the profile. It can be found in `name_odb10/scores_cutoff`. And yes, the higher the better.
> 3 votes

---
Tags: busco
---


# Question
Title: REMC SNP annotations for transcriptome prediction - epiXcan
I have been trying to use a pipeline for predicting gene expression for a target gene (Although you can do it for many) - the name of this method is epiXcan. Here is a link to the paper:

https://www.nature.com/articles/s41467-019-11874-7#Sec22

One of the very first steps in the pipeline is obtaining SNP 'annotations' from an epigenome reference consortium known as the REMC. In the methods paper they say:

"We first prepare eQTL statistics (computed with MatrixEQTL4) and SNP annotations (extracted from REMC https://egg2.wustl.edu/roadmap/web\_portal/). For each eQTL tissue, we use the matched REMC tissue to extract the corresponding annotations".

I have completed the first step of obtaining eQTL statistics. However obtaining SNP annotations is where I am lost.

In my case, I am only trying to predict the expression of one gene (DRD4) and so my region of interest is the pre-frontal cortex. If you visit the website they provide though, it doesn't seem as though annotating a large list of SNPs - at least in the way they are presenting - is possible.

Can anyone provide some insight or direction on how to proceed? How does one obtain SNP annotations from the REMC? I have tried emailing the authors of the paper and they have not been responsive.

Thank you,

# Answer
> 0 votes
After much research I found out how to obtain the annotations.

You must go to the REMC website: https://egg2.wustl.edu/roadmap/web\_portal/

Go to the chromatin-state-learning tab, and choose the n-chromatin state model desired - the core has 15 states, then 18, then 25 (imputed states).

Following this, you must go to the download URL link and, you will find a lot of files, primarily with the following coding information: E0\*.. This corresponds to a particular tissue/cell type and so you must choose the one relevant to your research question. The way they are correspond to each tissue/cell type can be found in the following paper:

https://www.nature.com/articles/s41467-019-11874-7#Sec22

---
Tags: sequence-annotation, snp, transcriptome, imputation
---


# Question
Title: Cannot blast database created with blastdb_aliastool: BLAST Database error: No alias or index file found
I want to merge/concatenate BLAST databases as suggested here.

This is my example script that explains what I'm trying to achieve better than words can:

```
#!/bin/bash

for file in FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn
do
    echo "## Make DB for $file"
 makeblastdb -in $file -dbtype nucl &>/dev/null
    blastn -db $file  -query query_single_nucl.fasta -outfmt 6
    echo
done

echo "## Create testdb.nal..."
blastdb_aliastool -dblist "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -dbtype nucl -title test -out testdb

echo
echo "## Test testdb.nal..."
blastn -db testdb.nal -query query_single_nucl.fasta -outfmt 6

```

This is the output:

```
## Make DB for FAM1079.ffn
test_query  FAM1079_00001   100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM3228.ffn
test_query  FAM3228_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM6161.ffn
test_query  FAM6161_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM19036.ffn
test_query  FAM19036p_000001    100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Create testdb.nal...
Created nucleotide BLAST (alias) database testdb with 12692 sequences

## Test testdb.nal...
BLAST Database error: No alias or index file found for nucleotide database [testdb.nal] in search path [/path/to/files::]

```

While the individual databases work, the concatenated databases don't.

# Answer
The problem was that the database must be specified without suffix, i.e. `blastn -db testdb` as opposed to `blastn -db testdb.nal`.

By the way, it produces exactly the same output as

```
blastn -db "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -query query_single_nucl.fasta -outfmt 6

```

and the performance is not better.
> 2 votes

---
Tags: blast, makeblastdb
---


# Question
Title: Tools for modelling and visualizing growth of cells
I am looking for a program/GUI app/package that would help me do simulations for cell growth. I have a microscopic mathematical model, and the scenario is basically the following: I start with a rod-shaped cell characterized my its center and direction angle, it grows until it reaches a certain length, it divides into two identical cells whose angles change randomly and the process continues.

Is there any program or package that I can use for such a situation? I have found several on GitHub, but they were mostly for very complicated biological processes. I use Python/Julia/Matlab and can code, but I was wondering whether there are any tools for such simpler scenarios.

# Answer
> 1 votes
If the scenario is simpler, and you can't find anything by searching, I would consider just writing your own.

The modelling can be done in any language but if compute time is an issue - for example you want to run a simulation with many cells or many simulations - then I would use a language that can support that such as Julia or NumPy in Python. In fact Julia is a nice language for this type of thing as you end up with readable code that runs fast.

For visualisation in Julia specifically you could use the animation functionality of Plots.jl, which is basic but there, or could try Makie.jl. See the Molly.jl documentation for some examples of visualising particle-type processes using Plots.jl.

---
Tags: python, modelling
---


# Question
Title: NCBI SRA database sample: control vs test
I was trying to download some data from NCBI SRA **(SRA059451)**. There are 27 samples available for **SRA059451**. But i am unable to understand which samples are **'control'** and **'test'** samples. Please help me regarding how to differentiate the samples.

# Answer
> 3 votes
You can download the information using this link:

I did it in R:

```
tab = read.delim("https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=SRA059451&result=read_run&fields=run_accession,experiment_title,sra_ftp,sample_alias")

head(tab)
  run_accession
1     SRR585697
2     SRR585700
3     SRR585726
4     SRR585727
5     SRR585728
6     SRR585729
                                                                experiment_title
1 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 171
2 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 238
3 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 241
4 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 242
5 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 243
6 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 245
                                      sra_ftp
1 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585697
2 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585700
3 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585726
4 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585727
5 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585728
6 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585729
                                    sample_alias
1     Symptomatic atherosclerosis gut metagenome
2 Symptomatic atherosclerosis gut metagenome 238
3 Symptomatic atherosclerosis gut metagenome 241
4 Symptomatic atherosclerosis gut metagenome 242
5 Symptomatic atherosclerosis gut metagenome 243
6 Symptomatic atherosclerosis gut metagenome 245

```

Normally the information will be here, but I cannot find what does metagenome 238 etc refers to. No table in the paper that can link sample to phenotype..

Most likely you need to contact the authors. This list might be useful though.

---
Tags: sequence-analysis, phylogenetics, metagenome
---


# Question
Title: MaCconkey Agar colour change in non-lactose fermentors
On Macconkey agar, lactose fermenting bacteria turn the agar bright red. Why does non-lactose fermenting bacteria turn the MaCconkey agar yellow?. Does it have to do with the alkalinity of waste products?

# Answer
McConkey agar is already red. It's trick is that it has two carbon sources. The bacteria that cannot use the more energetically useful lactose use the amino acids (peptone) generating ammonia which changes the pH and the crystal violet goes yellow. The ones that can use lactose acidify the agar a bit so go a bit more red —but not much.

This being the SE for bioinformatics, it should be said that if the genome of the organism in question is sequenced and deposited in NCBI you can find if the organism possesses the enzyme that allows the hydrolysis of lactose (beta-galactosidase) by doing a protein BLAST. The lacZ gene encoded protein in E. coli is NP\_414878 in NCBI. This goes in the accession number input box. The species in question goes in the organism box. If you get something with good coverage (80%) with more than 20% homology the organism is likely lac+.
> 1 votes

---
Tags: microbiology
---


# Question
Title: How to do BLASTP for short sequences (<20 aa) effectively?
I am trying to find sequence homology between viral sequences and my protein of interest. I have the sequences of their epitopes which varies from 5 to 500 amino acids long. For shorter sequences, it is recommended to use `blastp-short` but I am confused until which length it is okay to run `blastp`?

For shorter sequences, what are the parameters needs to be added for effective results?

`blastp -task blastp-short -query POI.fasta -db virus.fasta -out output.txt -outfmt 6 -evalue 20000`

# Answer
> 1 votes
You should increase the `-evalue 10` default option to a higher value because this value it's adjusted to the sequence database size. Check this table, maybe this can help you to decide which parameters are the best option for you.

```
    OPTION      TASK          DEFAUL    DESCRIPTION
    word_size   blastp        3         Word size of initial match.
    word_size   blastp-short  2         Word size of initial match.
    gapopen     blastp        1         Cost to open a gap.
    gapopen     blastp-short  9         Cost to open a gap.
    gapextend   blastp        1         Cost to extend a gap.
    gapextend   blastp-short  1         Cost to extend a gap.
    matrix      blastp        BLOSUM62  Scoring matrix name.
    matrix      blastp-short  PAM30     Scoring matrix name.
    threshold   blastp        11        Minimum score to add a word to the BLAST lookup table.
    threshold   blastp-short  16        Minimum score to add a word to the BLAST lookup table.
    window_size blastp        40        Multiple hits window size, use 0 to specify 1-hit algorithm.
    window_size blastp-short  15        Multiple hits window size, use 0 to specify 1-hit algorithm.

```

---
Tags: sequence-alignment, blast, covid-19, sequence-homology, blastp
---


# Question
Title: Possibility to save output blastn table in memory using biopython
Is there a possibility, using biopython to save output table of blastn in memory, not in file on hard drive, for process it using pandas for example and then delete it from memory?

I mean, is there a possibility to redirect output of "out" argument of this function into memory? In pandas DataFrame for example.

```
NcbiblastnCommandline(query="fasta1.fasta", subject="fasta2.fasta", out="output_table.csv", outfmt=6)

```

I always saved output tables in file before, but now I need to align one set of small fasta files against another set of small fasta files and I think that the best way is avoid writing of many small tables in hard drive because after recieving one table I have to read and process this table immediately and after that I don't need this table anymore.

# Answer
> 6 votes
You can specify send to stdout using `out='-'` in the Biopython wrapper.

```
from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd

cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()

rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', 'length',
        'mismatch', 'gapopen', 'qstart', 'qend',
        'sstart', 'send', 'evalue', 'bitscore']

data_types = {'pident': float, 'length': int, 'mismatch': int,
              'gapopen': int, 'qstart': int, 'qend': int,
              'sstart': int, 'send': int, 'evalue': float,
              'bitscore': float}

df = pd.DataFrame(rows, columns=cols).astype(data_types)

```

---
Tags: python, blast, pandas
---


title: SRA run data access

text:
# Question
Title: SRA run data access
Trying to download run data of Run accession **SRR2155174**. Here's the link: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR2155174 Please go to the link and go to the "Data Access" tab. Hereunder header "Name" you can find two links: https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR2155174/SRR2155174.1 https://sra-downloadb.st-va.ncbi.nlm.nih.gov/sos2/sra-pub-run-6/SRR2155174/SRR2155174.1 **Do these two links download the same data or they are for different purposes?**

Thanks

# Answer
It seems just like two different mirrors contains the same SRA archive. If you still doubt, you can download both (its have relatively small size) and compare their md5sum.
> 2 votes

---
Tags: sequence-analysis, phylogenetics
---


# Question
Title: How to simulate phenotype from real genetic data for GWAS purpose?
I'm trying to simulate binary phenotypes from the 1000 Genome Phase 3 datasets using gcta64 --simu-cc, but no success. Everything seems to be going well, but in the end I get:

```
Error: can not open the file [] to read.
An error occurs, please check the options or data

```

And the log shows:

```
Accepted options:

```

Here's the commands I'm using:

```
# Convert the VCF to plink format:
$ ./plink2 --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz --make-bed --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

# Try to simulate the phenotype:
$ ./gcta64 --bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs --simu-cc 500 500 --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

```

Here's the whole steps with the outputs:

```
$ ls
gcta64  plink2

$ ./plink2 --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz --make-bed --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
PLINK v2.00a2.3 64-bit (24 Jan 2020)           www.cog-genomics.org/plink/2.0/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log.
Options in effect:
  --make-bed
  --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
  --vcf ../../ALL.phase3.biallelic-only.vcf.gz.10kSNPs.vcf.gz

Start time: Sun Mar 22 11:56:23 2020
16384 MiB RAM detected; reserving 8192 MiB for main workspace.
Using up to 8 compute threads.
--vcf: 220000 variants scanned.
--vcf: ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pgen +
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pvar +
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.psam written.
2504 samples (0 females, 0 males, 2504 ambiguous; 2504 founders) loaded from
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.psam.
220000 variants loaded from
ALL.phase3.biallelic-only.vcf.gz.10kSNPs-temporary.pvar.
Note: No phenotype data present.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam ... done.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim ... done.
Writing ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed ... done.
End time: Sun Mar 22 11:56:28 2020

$ ls
ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed    ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam    gcta64
ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim    ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log    plink2

$ ./gcta64 --bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs --simu-cc 500 500 --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out ALL.phase3.biallelic-only.vcf.gz.10kSNPs
*******************************************************************
* Genome-wide Complex Trait Analysis (GCTA)
* version 1.93.0 beta Mac
* (C) 2010-2019, The University of Queensland
* Please report bugs to Jian Yang <jian.yang@uq.edu.au>
*******************************************************************
Analysis started at 11:59:16 GMT on Sun Mar 22 2020.
Hostname: Brunos-MBP

Accepted options:
--bfile ALL.phase3.biallelic-only.vcf.gz.10kSNPs
--simu-cc 500 500
--simu-hsq 0.5
--simu-k 0.1
--simu-rep 3
--out ALL.phase3.biallelic-only.vcf.gz.10kSNPs

Reading PLINK FAM file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam].
2504 individuals to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.fam].
Reading PLINK BIM file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim].
220000 SNPs to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bim].
Warning: Duplicated SNP ID "rs145607083" has been changed to "rs145607083_5264"
.Warning: Duplicated SNP ID "rs145607083" has been changed to "rs145607083_5265"
.Warning: Duplicated SNP ID "rs71955229" has been changed to "rs71955229_27061"
.Warning: Duplicated SNP ID "rs71955229" has been changed to "rs71955229_27062"
.Warning: Duplicated SNP ID "rs71589472" has been changed to "rs71589472_42505"
.Warning: Duplicated SNP ID "rs563156514" has been changed to "rs563156514_49111"
.Warning: Duplicated SNP ID "rs563156514" has been changed to "rs563156514_49112"
.Warning: Duplicated SNP ID "rs539504239" has been changed to "rs539504239_79196"
.Warning: Duplicated SNP ID "rs35739849" has been changed to "rs35739849_105514"
.Warning: Duplicated SNP ID "rs148795567" has been changed to "rs148795567_123134"
.Warning: Duplicated SNP ID "rs143101359" has been changed to "rs143101359_201815"
.Reading PLINK BED file from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed] in SNP-major format ...
Genotype data for 2504 individuals and 220000 SNPs to be included from [ALL.phase3.biallelic-only.vcf.gz.10kSNPs.bed].
Simulation parameters:
Number of simulation replicate(s) = 3 (Default = 1)
Heritability of liability = 0.5 (Default = 0.1)
Disease prevalence = 0.1 (Default = 0.1)
Number of cases = 500
Number of controls = 500

Error: can not open the file [] to read.
An error occurs, please check the options or data

$ head ALL.phase3.biallelic-only.vcf.gz.10kSNPs.log
*******************************************************************
* Genome-wide Complex Trait Analysis (GCTA)
* version 1.93.0 beta Mac
* (C) 2010-2019, The University of Queensland
* Please report bugs to Jian Yang <jian.yang@uq.edu.au>
*******************************************************************
Analysis started at 11:59:16 GMT on Sun Mar 22 2020.
Hostname: Brunos-MBP

Accepted options:

```

I'm opened to use other tools if you have any recommendation. Not sure, but looks like Plink doesn't do that (you can simulate phenotype as long as you also simulate the genetic data as well...)

* PS 1: I managed to simulate it before as answered here, but I've just realised it doesn't generate reliable ratios (although the logic in place, phenotypes are still too aleatory assigned - too low precision, recall and f1 scores).
* PS 2: Also asked in the BioStars forum

# Answer
> 1 votes
The question has been answered in the bioStars forum: https://www.biostars.org/p/428314/

It was missing the param `--simu-causal-loci`.

I created the `causal.snplist` file by finding out which SNP's are known causals of the treat I want to simulate. In my case, I want to simulate a Diabetes phenotype, thus I:

1. Got the known Diabetes causal SNP's from snpedia
2. Filtered out whichcausal SNP's also exist in my genomic data:

```
$ cat chr3_6_10.biallelic-only.bim | grep "rs4655595\|rs6718526\|rs4402960...

```
3. Created the `causal.snplist` with the causal SNP's as per documentation

```
$ head causal.snplist
rs1111875
rs12255372
rs1801282
rs237025
...

```
4. Run gcta-simu again:

```
$ ./gcta64 --bfile chr3_6_10.biallelic-only --simu-cc 500 500 --simu-causal-loci causal.snplist --simu-hsq 0.5 --simu-k 0.1 --simu-rep 3 --out chr3_6_10.biallelic-only.gcta-simu

```
5. Got the expected results:

```
$ head chr3_6_10.biallelic-only.gcta-simu.phen
0 HG00096 1 1 1
0 HG00100 1 1 1
0 HG00101 2 1 2
0 HG00102 1 1 2
0 HG00103 1 1 1

```

> The last 3 columns are phenotypes. 3 phenotypes were generated due the param `--simu-rep 3`.

---
Tags: plink, simulated-data
---


title: How to compute Shanon entropy

text:
# Question
Title: How to compute Shanon entropy
I want to analyze some viral sequences. I want to use shanon entropy as a measure to analyze the sequences. However I do not understand the concept of shannon entropy. Would anyone explain please how it works or advise on some resources i can use to understand shannon entropy.

# Answer
I think there needs to be clarification of the question. Shannon entropy is not the same thing as the Shannon-Wiener diversity index (H), commonly referred to as the Shannon index.

The reason I brought it up is because you have the tag "population genetics". Shannon diversity index is commonly used in population genetics to compare species richness (the number of species in a given area) and their relative abundance.

Shannon entropy is more to do protein structure, which isn't really population genetics, but can relate to how stable a given mutation might be. That is not really population genetics however, population genetics is about abundance of a given allele in this context.

I am pretty sure that on this site most responders will recommend alternative free energy approaches to perform this style of analysis.
> 2 votes

---
Tags: phylogenetics, phylogeny, statistics, genomics, population-genetics
---


title: Per Base Sequence Content in fastqc

text:
# Question
Title: Per Base Sequence Content in fastqc
I have a question regarding "Per Base Sequence Content" plot for "fastqc":

In the fastqc documentation, it is written: "In a random library you would expect that there would be little to no difference between the different bases of a sequence run, so the lines in this plot should run parallel with each other."

But I don't understand why different bases in a read should follow the same pattern of allele frequency ("A/T/C/G"). I mean they are different positions in the genome and it is normal that each position has different allele.

I would appreciate it if someone could help me, please.

# Answer
The positions you see in the x-axis of this plot are positions on your *reads*, not on the reference genome. Since each read comes from a random position on the genome, the frequency of A/T/C/G on each base position should reflect the base composition of the entire genome. Therefore, the percentage of each base should be approximately constant across the read.

Hope that helps!
> 2 votes

---
Tags: ngs, genome, fastqc
---


# Question
Title: what is the concept of N-channel geneChip?
I can understand two-color and one-color(Channel) gene chip; but there is another type of gene-Chip called "N-channel chips" ; I do not know how can several (\> 2) samples be hybridize in one chip. I cannot find some reference for the a material for this concept. See the following code that it has two channels "E“ and "Eb"; what does that channel mean? The expression matrices seem to have the same dimension. In my understanding, the two-color or two-channel microarray is to combine two samples (treated vs. control) in a chip, and the log(ratio) is used to express the differential expression.

```
AE1548 = getAE("E-MTAB-1548",
               path="/Users/zhang/Documents/2020/GEOsepsis/Data",
               local = TRUE, sourcedir = "/Users/zhang/Documents/2020/GEOsepsis/Data",
               type = 'full')
AE1548raw = ae2bioc(mageFiles = AE1548)#established from raw data; multi-channel data
AE1548expr1 = assayDataElement(AE1548raw,"E")
AE1548expr2 = assayDataElement(AE1548raw,"Eb")
  AE1548expr1[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1    273812.0     97733.5      205330
2        41.0        26.5          50
3        42.5        25.0          48
4        40.0        23.5          49
5        38.5        24.0          53
  AE1548expr2[1:5,1:3]
  SIRS_18exp2 SIRS_19exp2 SIRS_23exp1
1          25        18.0          36
2          25        18.0          37
3          25        18.0          38
4          25        19.0          37
5          24        18.5          37

```

# Answer
> 0 votes
**This is a mere long comment, which may help, but I do not know the details required for the answer.**

## Several (\> 2) samples hybridized in one chip

A traditional microarray has two parts, the mobile target, which is the labelled sample, and fixed to the filter paper, glass or doped silicon is the probe, which is not labelled and is in vast excess. You can label the target fluorescently with one of a variety of fluorophores, or with different radioisotopes, which have different energies. Nobody uses radioisotope microscopy for this except for the most obscure method, say nuclear run off assays, but I mention this just to cover all bases.

Whereas you can buy whatever fluorophore you wish, the detecting machine needs to be set up for it with lasers, beam splitters and filters. Most florescent confocal microscopes use a LED lamp, a splitter and a loads of cubes of filters. In the case of flow cytometry, you can get over 16 fluorescent channels at the same time. However, deconvolution becomes necessary when one channel "bleeds" into an other —I am not sure how one deals with it with a microarray, I suspect before the step to remove regional biases (the step with LOESS regression or similar).

## Affymetrix Gene chip (R)

Now, in the case of an Affymetrix Gene chip, I do not know if there is a multichannel machine, but it could be custom build. The spots on an Affymetrix chip are 3-5 µm, so regular microscopy optics still work. Although the silicon base is not transparent like a home made slide or an Agilent microarray (if these are still made). Hopefully, someone else can answer about Affymetrix machinery!

## Statistics

I am not familiar with why or how you would analyse a 3+ channel experiment. One guess could be to do a dose-response experiment. But there will be some problems:

* The dye flip correction would be an issue for greater than 2 samples
* Channel deconvolution may be required for dyes with overlapping spectra.

---
Tags: gene, chip-seq
---


# Question
Title: How do I differentiate outliers from in-group variations in a DESeq2 PCA plot of 9 samples distributed into 3 conditions?
I have a PCA plot from DESeq2's `plotPCA(vsd, intgroup=c("conditions"))` function. I have 9 samples distributed in to 3 groups of 3 biological replicates each. My reason for trying to visualise all on the same PCA plot is to see that samples from the same group (condition) cluster together and to see a separation between the groups.

The PCA plot I have is does not appear as I expect:

1-Can I just remove a sample which I deem is the outlier and perform a reclustering of the rest on the PCA?

2-What parameter can I use to declare that a point is an outlier rather than an in-group Biological variation?

3-I also wish to find out if it is allowed to plot more than two groups of sample at a time on a PCA plot.

Thank you in advance for your kind hints.

# Answer
> 1 votes
Since you have so few points here, it is hard to tell which are outliers by just looking at the PCA plot. I suggest you use correlation heatmap for your gene expression and check whether there are obvious outliers.

As for using more than one groups when plotting PCA, I think you can just try passing a vector with many groups to the argument `intgroup`.

---
Tags: rna-seq, deseq2
---


# Question
Title: What are the largest fully publicly available molecular QTL datasets?
I'm familiar with the sQTL and eQTL GTEx data that can be downloaded from the GTEx Portal, but I'm interested in gathering at least two other types of molecular QTL datasets:

1. Larger eQTL and/or sQTL datasets (greater overall sample size, even if performed only for blood cells)
2. Other functional molecular QTLs such as methylation QTLs, histone QTLs, DNA accessibility QTLs, etc.

# Answer
> 1 votes
The largest eQTL datasets I can think of is probably eQTLGen. This is a meta-analysis of many studies using \>30,000 samples from whole blood. This might be less ideal if you want cell-type-specific QTL information.

For other molecular QTLs, I do not think there are large-scale datasets yet.

---
Tags: genomics, public-dataset, qtls
---


# Question
Title: How to identify latent variables in single-cell RNA-Seq data
I have a single-cell RNASeq sample, in which I'd like to identify latent variables (e.g. response to stress) that I think might be affecting the clustering.

The approach I was planning to use is to reduce the dimensionality of the dataset using one algorithm among:

* PCA
* ICA
* NMF

Then identify the most important biological processes for every component (e.g. in each Principal Component), and then finally regress out the processes that might not be relevant to the analysis or that might affect the downstream clustering.

I have a few questions:

* Are there issues in this approach?
* Is there a collection of gene sets (similar to MSigDB) that can be used to identify noise latent variables?
* Which algorithm (PCA, ICA, NMF, or others) is more suited to this kind of analysis?

# Answer
> 1 votes
I think to identify latent variables, PCA is probably not going to work. NMF might be worth trying. You might want to check out a method called consensus NMF (cNMF) (https://elifesciences.org/articles/43803) published in eLife. The author described some latent variables corresponds to cell types, while others corresponds to cell-cycle or metabolic states of cells.

---
Tags: scrnaseq, clustering, pca, gene-ontology, noise
---


title: Joinning columns from different files

text:
# Question
Title: Joinning columns from different files
I have 48 files in which I have

r1

```
Organ-specific cancer signature Contribution
Esophagus_A 2948.114
Esophagus_B 815.6483
Esophagus_C 1278.381
Esophagus_D 883.5358
Esophagus_E 1446.103
Esophagus_F 880.0416
Esophagus_G 1262.907
Esophagus_H 743.2682

```

r2 .. r15

and

Non\_res1

```
Organ-specific cancer signature Contribution
Esophagus_A 2660.722
Esophagus_B 1471.017
Esophagus_C 1639.904
Esophagus_D 0
Esophagus_E 2580.775
Esophagus_F 946.3675
Esophagus_G 2013.654
Esophagus_H 0

```

...

Non\_res33

I want a unified file in column one is Organ-specific cancer and 48 more columns are the second column of each separate file

How I can do that?

# Answer
> 1 votes
If the first column is always the same:

If the first columns in not same:

---
Tags: r, bash, awk
---


title: Warning in fastqc

text:
# Question
Title: Warning in fastqc
I am checking the quality control of my sequences using the Fastqc tool. For some steps, like "Per base GC content", etc., I received a warning. So, I was wondering whether I Should also take care of warning messages in the Fastqc tool or I should only care about errors! I would appreciate it if someone could help me with this, please.

# Answer
> 1 votes
Generally yes, but it depends the type of study. Remember that Fastqc was developed for bacterial genomes, mainly, so an unusually shaped distribution in that particular dataset could indicate a contaminated library. But, for example, in an amplicon 16S amplicon dataset it's expected that other modules fail (Duplicate Sequences, Overrepresented Sequences, etc.) not only GC.

---
Tags: ngs, genome, sequence-analysis, genome-sequencing, fastqc
---