# Question
Title: Finding common and unique data set by comparing two files based on their column and to split the columns multiple strings to print in output
I have very large sizes tab-delimited .vcf files and want to match these two / or 3 files based on their position and print to a new .csv file

**File structures:**

File\_1: tab-delimited file (.vcf) and its as column names as follows `(line number 3439) #CHROM POS ID REF ALT QUAL FILTER
 INFO`

**File 2: same as file\_1 column names**

`(line number 3407) #CHROM POS ID REF ALT QUAL FILTER
 INFO FORMAT SAMPLE_1 SAMPLE_2 .....`

In file\_2, column 7 (INFO) contain many substrings like `AC=46,2;AF=0.958,0.042;AN=48;DP=269;ExcessHet=3.0103;FS=0.000;InbreedingCoeff=0.5411;MLEAC=92,4;MLEAF=1.00,0.083;QD=25.36;SOR=2.488 from these strings , have to print only information of AC=, AF=, AN=, DP=`

**Desired output files to generate:**

1. position matched in two of files: common\_position\_matched.csv

if I have more than 3 files also, the output should be one file and important thing, if a position column (1-POS) only matched only in 2 files and 1st file line should be NA NA NA

```
file_1  CHROM  POS    REF  ALT     file_2 #CHROM  POS      REF  ALT  INFO

1       22    10511521  C     T         1    22  10511521   C    T   AC=46,2 AF=0.958,0.042 AN=48 DP=269
2       22   10510544   G     A         2    22  10510544   G    A   AC=49,2 AF=0.958,0.042 AN=89 DP=536
3       22   10515068  AGAT,T AGAT,AT   3    22  10515068   AAA  AAAGG,A,GAA AC=100  AF=0.958,0.042 AN=62 DP=123
4       22   10515118  A G,   TAA       4    22  10515118   AG,  TAA AC=32   AF=0.958,0.042 AN=45 DP=500
5       22   10515118  AAAG   A         5    22  10515118   AATG A   AC=50   AF=0.958,0.042 AN=49 DP=129

```

**note:** while doing matching, not removing the duplicates, because in the same position there may be an addition or sometimes it may be deletion.

2. unique position of each file, in tab-delimited

output: File1\_unique.csv and File2\_unique.csv etc.

so far was able to read the file and match them according to position and print the output, but was not able to write efficient code

import pandas as pd
df1 = pd.read_csv("File1_3.vcf",sep='\t',usecols = ['POS']) ## Reading file1
df2 = pd.read_csv("file2_3.vcf", sep="\t", usecols = ['POS']) ## Reading file2
df3 = pd.concat([df1,df2], sort=True) ## Combining both the dataframes
df4 = df3.drop_duplicates(keep=False) ## Dropping the duplicates (intersect)
df4.to_csv("c3-UniquePosition_of_bothData.csv", sep="\t", index=False, header=True) ## Writing the unique to both
df1_Uni_file1_c3Posi = pd.merge(df4, df1, on='POS', how='inner') ## Identifying the unique position of File1
df2_Uni_File2_c3Posi = pd.merge(df4, df2, on='POS', how='inner') ## Identifying the unique position of File2
df_File1_File2_common_c3Posi = pd.merge(df1, df2, on='POS', how='inner') # Identifying the common chr-position of File1 and File2```

Program 2: (giving original file without editting)
import pandas as pd
df1= pd.read_csv("File1_22.vcf.gz", sep="\t", skiprows=3438, usecols = [0,1,2,3,4])
df2 = pd.read_csv("File2_22.vcf.gz", sep="\t", skiprows=3406, usecols = [0,1,3,4,7])
#writing the output files
#df1.to_csv("File1_c22.csv")
#df2.to_csv("File2_c22.csv")
#mergeing
df3 = pd.merge(df1, df2, on='POS', how='inner', sort=True)
df3.to_csv("common_position.csv", sep=",", index=False, header=True)
#df3 = pd.concat([df1,df2], axis=1).to_csv('check1.csv') # this command join multiple output to single output

Could any one give efficient python pandas script, to do this

Thanks All

# Answer
> 0 votes
One way to do this is to use a dictionary. You can read the first file line by line, and for each line, use "chr:position" as a key, and the corresponding line as the value. If you want to handle duplicates, you can make the value a list and store multiple lines with duplicated "chr:position" as elements in that list.

After constructing this dictionary, read the second file line by line. For each line, set "chr:position" as the so-called name of that line. Then using this "chr:position", you can easily find the corresponding line in the dictionary. Append the line from the second file to that line from the first file.

After finishing reading the second file, you are ready to write the new file to output. To do so, iterate over the keys in the dictionary, and write the value to each key to a new file. Remember that each value in that dictionary is now a new line.

I choose to use a dictionary because they are basically a hash table, meaning that using the key you can instantly read that value in the dictionary, instead of iterating every element in a list to find that element. This is very handy when you have a large dictionary.

---
Tags: python, perl
---


# Question
Title: Is there a computational tool or possibility to identify mRNA isoforms from the count matrix of a bulk RNA sequencing dataset?
I have the counts matrix of an RNA sequencing dataset of fibroblasts and I wish to identify isoforms of a particular gene of interest in it. Can anyone please hint me on a bioinformatics method to achieve this?

In literature, 9 different isoforms of my gene of interest have been reported. My colleague who does the wet lab experiments has identified 3 of those isoforms through western blotting in the fibroblasts that later got sequenced. There might be more isoforms in these fibroblasts.

I will be grateful for any computational hint to identify gene isoforms from my dataset, irrespective if they are differentially expressed or not.

Thank you in advance for your kind help.

# Answer
You can back to your bam file and use cuffdiff tools it is very able to do your job follow this to liks http://cole-trapnell-lab.github.io/cufflinks/manual/

https://github.com/Jeanielmj/bioinformatics-workshop/wiki/Differential-Analysis-with-Cuffdiff

My best
> 2 votes

# Answer
If literally all you have is the count info by gene, you can't magic up the isoform breakdown.
> 1 votes

---
Tags: rna-seq, deseq2, rna, rna-alignment, rna-splicing
---


# Question
Title: Overlap in Paired-end Reads for Sequencing?
Regarding Sequencing:

Do/can the paired-end reads have overlaps?

# Answer
Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.
> 4 votes

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


title: Location of Reads in Sequencing

text:
# Question
Title: Location of Reads in Sequencing
I have questions regarding sequencing:

Are paired-end reads from different stands? What about single reads, are they coming from both strands?

# Answer
> 2 votes
Mates in a pair are giving signal for a single strand, though this is only apparent in bisulfite-treated samples. The reason for this is that you are sequencing the ends of fragments of material and it's random which strand of this (or both in different clusters) will get sequenced. Different pairs of reads will then randomly sequence one strand or the other.

Single-end reads are the same as paired-end reads, the second read in the pair is just not made.

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


# Question
Title: GWAS, MWAS, EWAS: what are the (in)dependent variables?
I started reading some papers on X-wide association studies, where X can be metabolome, epigenome, etc... The authors usually describe which are the dependent and which are the independent variables of their models. Taking as an example a metabolome-wide association study (MWAS), the authors use as model the following: `metabolite` ~ `condition`. That is, each metabolite is considered the dependent variable and the condition (e.g., smoker vs. non-smoker) is considered the independent variable. I was a little bit confused, so I looked into an example study using GWAS: in this case the model can be written as `condition` ~ `SNPs`, which is what I would expect (knowing nothing about GWAS but something about Data Science) and quite the opposite from the MWAS described above.

Does the choice of the dependent variable depend on the questions asked in the study, or does it exist some sort of *agreement*? Moreover, regarding the GWAS model described above, I read that usually researchers fit a model for each SNP: what's wrong with fitting a model which consists of several independent variables (i.e., `y ~ SNP_1 + SNP_2 + ... + SNP_N`)?

# Answer
> 3 votes
To answer the first part of your question, the dependent and independent variable of X-WAS is kind of arbitrary and dependent on the question you asked. But it gradually becomes a convention in the field after the initial name and concept are accepted by the community. For example, GWAS from the very beginning is written as `condition ~ SNP`, and there is no doubt about it. Later, people did something like `condition ~ gene expression` and called it TWAS, which is also widely accepted now. In your case, MWAS is indeed `metabolite ~ condition`, and the only important thing is to know the convention in the research community so people are on the same page when they talk to each other.

As for why running GWAS using single SNP. In a GWAS, it is typical to use 500K to 7M SNPs. So this is a large number. It is obvious that we cannot include all the SNP in one single model, then the question becomes how many SNPs to include in each model, and how to choose those SNPs. This is very hard to decide. Also, due to link disequilibrium, many SNPs are correlated, which means co-linearity, a thing you want to avoid in linear models. It is not clear how to correctly interpret the effect sizes in a model with many SNPs as well. Therefore, due to feasibility, simplicity, and interpretability, it might be best to perform regression on every single SNP.

Now since we are running GWAS in an SNP-by-SNP fashion, there are a lot of work that needs to be done after GWAS in order to find the true causal SNP. Many fine-mapping methods have been developed to tackle this problem.

Interestingly, there are some multi-variate GWASs methods, but they do not come in the flavor of simple linear models. Instead, sometimes complex Bayesian statistics are used, for instance.

---
Tags: gwas, linear-regression, logistic-regression
---


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

```
CONFIRMED_DEATHS/CONFIRMED_INFECTIONS

```

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

```
CONFIRMED_DEATHS/(CONFIRMED_DEATHS + CONFIRMED_RECOVERIES)

```

Which formula is applicable to which situations?

# Answer
This may not strike most as a bioinformatics, but getting the key clinical outcome is essential in understanding the molecular basis of pathogenicity.

I think the mortality rate is over-reported. This is not to say the situation of 2019-nCov is not serious - it is very serious.

The two essential factors missing in your equation for 2019-nCov are:

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

Old people have a weakened immune system (e.g. reduced T-cell response). The very young have an underdeveloped immune system. The problem is hospital cases are skewed towards those with severe infection, which is likely the old and very young. In 2019-nCov the first infecteds were middle aged men, but likely to be restricted to those exposed to the "Wuhan meat market".

**Problem** Asymptomatic "spreaders" are known in influenza and are suspected in 2019-nCov. Asymptomatics will inflate the mortality rate (because we're not detecting them).

**Solution(s)**

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

I would assume the hospital is the place virus-induced mortality is reasonably accurate.

The second formula you described works very well for Ebolavirus, no-one has ever heard of asymptomatic Ebolavirus. It would also work well for another betacoronavirus Middle Eastern respiratory syndrome. The key problem with 2019-nCov is we don't know the number of asymptomatic patients or patients with such mild infections they never report illness.

**Summary** In summary, I think the mortality is being over-reported if it is based on your equation 2. Equation 1, however, will also over-report the mortality rate.
> 4 votes

# Answer
To add to Michael G.'s answer, there is a further wrinkle in the case of 2019-nCov that is due to the exponential growth in confirmed infections daily.

That growth affects your second formula `CONFIRMED_DEATHS/(CONFIRMED_DEATHS + CONFIRMED_RECOVERIES` since deaths probably occur earlier on in the disease trajectory than recoveries do. Since the number of confirmed cases is growing exponentially on a daily basis, the number of deaths will presumably be inflated relative to the number of recoveries at this point.

This feature may have the opposite effect on your first formula `CONFIRMED DEATHS/CONFIRMED INFECTIONS`, depending on when the infections are confirmed. If the infections are generally confirmed several days before death occurs, then a large portion of confirmed infections are still too early to have ended in death and `CONFIRMED DEATHS/CONFIRMED INFECTIONS` actually underestimates the actual number of fatalities.
> 2 votes

# Answer
I was facing this problem in January when coming up with a methodology to forecast on the UrbanSurvival.com site.

My solution - using low data collection was to use a daily rolling rate of infections and divide that into the 5-day previous reported death counts.

This was based on the assumption that, for a given cohort, the people who get sick today will "get dead" in about 10-days time.
So half of this (5 days) ought to get me into the ballpark.

In Feb, when CDC was reporting 3.4% mort. rate, we'd been reporting that for a couple of weeks.

As a formula it's not unlike deriving an exponent.

That is: Cases Today (T1) / Deaths 5-Days Ago (T2)

Unfortunately, this is presently running a bit higher mortality rate now (4.5%) but I would expect on the backside of the expansion of cases it will decline.

Here's hopin'...
> 2 votes

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


# Question
Title: Feature extraction methods that can handle inconsistent numbers of atoms for molecular dynamics
I want to compare the protein dynamics ) pH 7 versus pH 3, or ) wild type versus mutant

The protein will have slightly different number of atoms at each condition, due to protonation or mutation, which is bad for machine learning. The alternative is to only consider the alpha-carbon atoms, but this will lose significant amount of information (e.g. side chain, salt bridge).

Is there a way in molecular dynamics to extract features from ensembles of different number of atoms?

Note: I have used common methods to decipher the difference between two conditions, including RMSD, RMSF, Rg, secondary structure, salt bridge, SASA, etc. They can explain the difference to some extent. However, I feel the selection of properties is "subjective", and there must be some hidden features not known. So I am seeking a machine learning method, that can "automatically" capture the most essential properties based on the coordinates as inputs.

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

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

## Heavy atoms are common

Assuming the pH difference only results in protonations, the protons differ between the two conditions, but not the heavy atoms. These will retain all the side chain information.

## Merit of backbone only

*(Section no longer relevant due to OP edit)*

The common metric used to compare two structures is root mean square deviation (RMSD) of the alpha carbons to assess how much does the backbone differ and with analysis of MD simulations it is used to assess whether one is more dynamic than the other (e.g. this Bioinformatics SE question). The reason why the RMSD is done with the alpha carbons and not the minimum common subset of atoms is that backbone changes are large scale, while the side chain changes only affect locally —so their inclusion would needlessly dilute the score. The C prime and backbone nitrogens behave like the C alpha so there is no need for these.

## Extraction of coordinates

Assuming by features you mean atomic coordinates, there are loads of ways. Using Pymol module for Python, one could do something like:

```
import pymol2
with pymol2.PyMOL() as pymol: #can be parallelised
    pymol.cmd.load('trajectory.gro', 'mytraj')
    pymol.cmd.remove('h.') #remove hydrogens
    pymol.cmd.intra_fit("name CA", 1) #align all steps in the trajectory
    for i in range(1, pymol.cmd.count_states('mytraj')):
        for atom in pymol.cmd.get_model(f'state {i}').atom:
            atom.coord #List[int int int]

```

One could change `pymol.cmd.get_model(f'state {i}')` to select `state {i} and name CA` for Cαs only.

## Edit: Mutated residues

In the case of mutated residues or residue differences due to the two proteins being close homologues, "mutating" the differing residues to alanine, i.e. by stripping all heavy atom beyond the Cβ (keeping C′, Cα, H, N, O & Cβ) in the MD trajectory —or glycine if one is glycine. If they are homologues with different lengths, this approach would mean skipping residues that match up with gaps in the MSA. Even though not uncommon for statistically analyses, it does pose some issues statistically. It may unbury hidden residues, but length variant generally happens with loops.

In the snippet above, the pymol mutagenesis wizard could be used, but that is overkill, does not work in parallel and the approach does not port to other PDB reading modules in whatever language one wants. So say you have a list `to_be_trimmed` of integers of PDB residue indices (not pose residue indices):

```
if atom.resi not in to_be_trimmed or atom.name.strip() in ('C', 'CA', 'H', 'N', 'O', 'CB'):
      do whatever with atom.coord

```

In acidic environments, you get some funky things happening like exposed aspartate-glycine residue pairs forming isoaspartate or exposed aspartate-proline resulting in a cleavage, etc. For these cases, which are nearly never modelled, ignoring these values is the best option.

## Edit: small comment on SASA

The investigated property list is very extensive and some reader may adopt them. So I thought I should mention that per residue SASA is generally converted to relative surface area, which is residue independent, by dividing it by the maximum SASA for that residue (values from Tien et al. 2013: `'A': 121.0, 'R': 265.0, 'N': 187.0, 'D': 187.0, 'C': 148.0, 'E': 214.0, 'Q': 214.0, 'G': 97.0, 'H': 216.0, 'I': 195.0, 'L': 191.0, 'K': 230.0, 'M': 203.0, 'F': 228.0, 'P': 154.0, 'S': 143.0, 'T': 163.0, 'W': 264.0, 'Y': 255.0, 'V': 165.0}`), where RSA \> 0.2 is an exposed residue.
> 1 votes

---
Tags: machine-learning, features
---


# Question
Title: Since every human has a different DNA (different combinations of C, G, A, T) what does it mean to have the genome done?
I'm confused about the difference between genome and DNA. Is it correct to say that the same type of bacteria has the same DNA? But my understanding is that it is not correct to say that the same type of human has the same DNA, since every human has a different DNA. What am I missing here?

# Answer
All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are also often fatal, so the random mutations we see of living adults are much more restrictive than all the random mutations that occur.)

To identify the human genome is to study many people's DNA and to label the parts that are shared between everyone, the parts with two or three variations, and the parts with even more variations. Even though every person has different DNA, we can still say they fit the pattern, just like every T-shirt may be unique but they all fit the T-shirt pattern and not the trousers pattern.
> 21 votes

# Answer
DNA is only the substance that makes up the genome.

DNA can come with many forms. For example, the fragments you get by PCR in the lab are DNA. A bacterium also has DNA molecules in it. You may also find DNA fragments in many other places, like in our blood, or at a crime scene.

Genome, however, is a very specific term. It means all the heritable substances in an individual or an organism. Most times, when we talk about the genome in eukaryotes, we refer to the nucleus genome, that is, the DNA in the nucleus.

Note that genomes do not have to be made of DNA. RNA viruses, for example, only have RNA as their genomes. Yet we still refer to that as a genome, because that is the heritable substances for those viruses.
> 8 votes

# Answer
DNA is the product you get if you paste the nucleotides A, T, C and G together. This is DNA and all genomes consist of these basic building blocks like a book consists of letters. Still, the genome, or the content of the book is the specific combination of these letters plus the proteins that bind to the chromosomes, the folding of the DNA and all the epigenetic modifications that can happen. So yes, we all have the same DNA in terms of the same building blocks, but the content is different.
> 4 votes

# Answer
Might have misunderstood you - but "we" have sampled and "decoded" the DNA of a few humans (often volunteers who surrendered their body to science after their death) - but this gives of course only the genetic code of those few individuals.

When scientists say we decoded the genome they mean they deciphered which parts of the DNA "do" specific things - the actual "genes" or gene-sequences - usually done when looking at the DNA molecule packed up into the chromosomes to give easier "orientation".

And by my understanding the first task was accomplished - supported by use of supercomputers, but the deciphering of the gene-sequences still has a long way to go for scientists to state "complete." But there is growing success - And as usual scientists will not stop until complete.
> 1 votes

# Answer
"To have the genome done" means that we have sequenced all the DNA in particular reference individuals. You are correct that every human has different DNA, so there is no such thing as a single "human genome", as (almost) everyone's is unique. But the vast majority of everyone's DNA is identical - that's what makes us human. Having the full sequence of a few "normal" individuals who don't have glaring genetic defects is extremely useful for learning about the genome of *all* individuals, since so much of it is the same from person to person.

I've used the term "reference" a few times, since we use a few individuals' genomes as an expectation of what a human genome looks like. The process of sequencing DNA involves breaking it into many bits, and then recombining them using a computer. You can put the pieces together without the reference, but it is much more difficult and time consuming. Even if there are differences from the reference, having an expectation of what \>99% of the genome should be makes the problem much simpler. In a way, it's like the difference between doing a jigsaw puzzle without knowing what it's supposed to look like, and using the picture on the box as a reference. Even if your exact jigsaw puzzle has a few different pieces, the box reference is mostly the same and immensely useful.
> 0 votes

# Answer
Think of it as a special book. Everyone's book is different, but they all have the same number of chapters and each chapter covers the same topic in each individual book (sorta). The problem is we can only read a few words at a time.

When the genome is "done", it means that someone has read a bunch of these sentence fragments from one person's book and then assembled all the words in the correct order to make a complete book. Then we can get a good idea of how many chapters it has and what those chapters are about.

However, the human genome wouldn't be a novel, it would be more like a set of encyclopedias.
> 0 votes

# Answer
When people say "the genome" for a species, they mean a reference genome. To quote wikipedia,

> A reference genome (also known as a reference assembly) is a digital nucleic acid sequence database, assembled by scientists as a representative example of a species' set of genes. As they are often assembled from the sequencing of DNA from a number of donors, reference genomes do not accurately represent the set of genes of any single person.

The genome is updated over time as more is learned, the most recent human reference is hg38.

For individualized results, the patient is sequenced and the data is compared to the reference genome. Then, the differences are presented as a list of genetic variants. This is more efficient than looking at the entire genome for each patient, because on average greater than 99.9% of sites match the reference.
> 0 votes

# Answer
Yes you are right. 99.9% of our genome remains same, it is only the 0.1% which makes each individual unique from the rest of us. Single nucleotide polymorphisms (SNPs) contributes to the 0.1% variation in our genome. SNPs explain the variation in eye color, susceptibility towards disease conditions etc
> 0 votes

---
Tags: genome, phylogenetics
---


# Question
Title: Are the conclusions in "The proximal origin of SARS-CoV-2" legit?
I don't have any background in genetics and bioinformatics, so I ask you if you think that the arguments provided in the article *The proximal origin of SARS-CoV-2* by Andersen *et al.* are convincing. In particular:

> *While the analyses above suggest that SARS-CoV-2 may bind human ACE2 with high affinity, computational analyses predict that the interaction is not ideal and that the RBD sequence is different from those shown in SARS-CoV to be optimal for receptor binding. Thus, the high-affinity binding of the SARS-CoV-2 spike protein to human ACE2 is most likely the result of natural selection...*

I interpret this passage as "if the virus had been engineered, they would have done a better job". Is this interpretation correct? And if it is, does it sound convincing to you as a proof that the virus developed in nature?

The second motivation given in the article is the following:

> *Furthermore, if genetic manipulation had been performed, one of the several reverse-genetic systems available for betacoronaviruses would probably have been used*

Again, do you think that this exclude the possibility of a human intervention in the creation of the virus?

Thank you for your opinions.

Edit: I am a fan of Occam razor, I know that the scenario in which the virus originated in the wild is by far more likely than the human engineered scenario. I just want to know if, given our current knowledge in genetics, would have been possible for some high skilled researchers to engineered COVID19 (without implying that they did this with any bad purpose)?

Edit 2: I share the concerns of the author of this post, that is, that the findings of the article are merely opinions and some of arguments are misleading. Anyone having a solid technical background can comment on this? *China owns Nature magazine’s ass – Debunking “The proximal origin of SARS-CoV-2” claiming COVID-19 definitely wasn’t from a lab*

# Answer
In summary, the authors are saying the complete opposite of "human intervention".

> While the analyses above suggest that SARS-CoV-2 may bind human ACE2 with high affinity, computational analyses predict that the interaction is not ideal and that the RBD sequence is different from those shown in SARS-CoV to be optimal for receptor binding.

The interaction with ACE2 is taking place on the Spike protein (S) in both COVID-19 and SARS-Cov. The S protein is very different between COVID-19 and SARS both in terms of sequence and COVID-19 includes an additional furin cleavage site, meaning structurally it is quite distinct. The "computational analysis predicts" is difficult because there is no crystal structure for COVID-19 S protein therefore modelling ligand binding is essentially guess work, particularly as it comprises an additional (furin) cleaved domain. What the authors appear to conclude is that they cannot map ACE2-Spike binding of COVID-19 using an established SARS S protein structure. I assume they are homology modelling SARS-CoV onto COVID-19 S protein sequence data. In my personal opinion there is no way to homology model an additional furin protease cleavage site (furin is a host protease) so it would be difficult for "good fit" to occur, but it is just my opinion and the authors may have established methods to overcome this notable limitation, e.g. similar phenomena occur in influenza.

> Thus, the high-affinity binding of the SARS-CoV-2 spike protein to human ACE2 is most likely the result of natural selection..

Their conclusion of "natural selection" is not strictly accurate and should be "Darwinian positive selection", but its a small issue. What they are saying is Darwinian adaptation has occurred as opposed to "purifying selection, i.e. conservation. Thus they are definitely excluding human intervention because that is not "natural selection"

> Furthermore, if genetic manipulation had been performed, one of the several reverse-genetic systems available for betacoronaviruses would probably have been used

They appear to be saying that if a synthetic virus was constructed it would be created using an established reverse genetics system, and all those produced to date look nothing like COVID-19.

The phylogenetics argument is always a stronger argument than either of the above (although my opinion isn't entirely objective). Phylogeneticists have had a long history in opposing conspiracy theories, e.g. polio vaccination, HIV and the list goes on, so it is much more natural territory. An example of such an argument is here which @terdon helpfully pointed out.

---

To address the questions of @Hans

**Question 1**

* Purifying selection = any evolutionary change is deleterious, so the virus is less able to transmit between humans. The vast majority of mutations are deleterious in classical Darwinian thinking\*
* Adaptation = the amino acid change occurs to ensure the virus is better able to transmit between humans. In phylogenetics this is called 'positive selection' which is detected at a nucleotide level, this is however a very stringest test.

\*, There is something called nearly neutral theory, but thats just complicated

---

**Question 2**

What we are talking about is a virus infectious clone system. For a 30kb virus this would not be trival at all, but I agree for 10kb viruses it much easier. However, the phylogenetics are opposed to this, because SARS-CoV-2 is 5% divergent from a virus (RaTG13) isolated a long time before the current epidemic. You can't just engineer a virus that is 5% divergent from a bat virus across the entire genome and never been seen before. It would be a work of genius the likes we have never seen using bioinformatics as yet unheard of, because what would happen is you would encounter endless deleterious mutations. Even if you got over all those, how do you know what fitness traits you are aiming for?

**Serial passage** This level of change is doable by serial passage of the virus, but this leads to attenuation of virus (less dangerous) not increased fitness. The yellow fever virus vaccine was created in this way. Even if a virus had been serially passaged we could tell, because we know the mutational patterns of other serially passaged viruses.
> 4 votes

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


# Question
Title: Analyzing proteins based on sequence similarity
Using NCBI Gene database, I retrieved list of 7015 protein which are implicated in autism pathogenesis. This list includes tax-ID, organism name, symbol, aliases, map location, chromosome, genomic nucleotide accession version, start position on genomic accession, end position on genomic accession, orientation, exon count and OMIM). Among these 7015, 1144 proteins are involved in neurodevelopmental functions. I want to know that how many of these 1144 proteins are functionally similar to 7015 proteins? Is there a method to query 1144 proteins against 7015 proteins and group them based on structural similarity/ function similarity/ gene family?. Since I have large set of proteins, I am trying to understand if sub-set of proteins share similar function with this large set of proteins.

Thanks

Priya

# Answer
I suggest you give DAVID a try. Specifically, their Functional Annotation tool. Just enter your list of protein IDs, and it will return groups of proteins where particular GO functions are overrepresented.

This should at least help you start categorizing your list into functional groups.
> 2 votes

---
Tags: structural-variation, gene-ontology, sequence-homology, multi-omics, paralogs
---


# Question
Title: How to decide number of neighbors and resolution for Louvain clustering
I am using Louvain clustering (1,2) to cluster cells in scRNAseq data, as implemented by scanpy.

One of the parameter required for this kind of clustering is the number of neighbors used to construct the neighborhood graph of cells (docs).

Larger values result in a more global view of the manifold, leading to lower number of clusters, while reducing the number of neighbors goes in the opposite direction. However, it is unclear how to choose this parameter.

The resolution parameter seems to work in the opposite way.

Do you know of any methodology and/or rule-of-thumb to define these parameters? E.g. depending on the size of the dataset?

1. Levine, Jacob H., et al. "Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis." Cell 162.1 (2015): 184-197.
2. Blondel, Vincent D., et al. "Fast unfolding of communities in large networks." Journal of statistical mechanics: theory and experiment 2008.10 (2008): P10008.

# Answer
> 1 votes
A general rule of thumb is that in order to improve the variance $n$ times you need $n^2$ neighbours. This is only applicable if you consider the $n^2$ nearest neighbours of a cell to be biologically identical (i.e. "similar enough"); if your data includes 10 types of cells with 10 cells each, then using the 20 nearest neighbours for smoothing will obscure the data.

As far as I know, there is no single best answer to this question. I would suggest trying different numbers and sticking to what agrees more with the biology of the dataset.

# Answer
> 1 votes
I read a paper where the authors used different combinations of parameters to iteratively perform LSI until some criteria. Maybe a similar strategy can be applied to Louvain clustering as well?

The paper can be found here: https://www.nature.com/articles/s41587-019-0332-7. You can also check their code to see how they implemented it.

---
Tags: scrnaseq, clustering
---


# Question
Title: How can I download the SNP genotyping information from GTEX?
I am new to the field of Bioinformatics. GTEx database is great; I know how to check whether a SNP is eQtl **with the online tools** of GTEx.

But my question is how to download the right data to draw a boxplot on my own. like this paper:

Please help me . I have been stuck on this for many days.

# Answer
There are two types of data you can get from GTEx. The first is summary statistics, basically the slopes and p-values you see in the boxplot. You can also get expression levels for all genes measured in GTEx.

However, in order to plot boxplot by yourself, you also need genotype data. This is not publicly available because of privacy issues. In order to get it, you have to apply to it on dbGaP. I think you can find relevant information on the GTEx website.
> 1 votes

---
Tags: snp
---


title: Merging based on Chromosome coordinates

text:
# Question
Title: Merging based on Chromosome coordinates
I have a simple task which seems very complicated now Two data-frame I have One contains the Chromosome coordinates of distal enhancer elements and other one rlog values of accessibility which I want to map.

The distal coordinates data

```
dput(head(a))
structure(list(Chr = c("chr6", "chr11", "chr20", "chr20", "chr10",
"chrX"), Start = c(154986726, 70314780, 15895696, 40706992, 29132014,
129525939), End = c(154987062, 70315059, 15895987, 40707264,
29132266, 129526265)), class = c("tbl_df", "tbl", "data.frame"
), row.names = c(NA, -6L))

```

The Rlog\_ATAC\_location data

```
Chr    Start    End  CMP1  CMP2  CMP3  CMP4  CMP5  CMP6  GMP1  GMP2  GMP3  GMP4  GMP5  GMP6  HSC1  HSC2  HSC3  HSC4  HSC5  HSC6
  <chr>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 chr1   16103  16348 2.70  2.99  2.87  1.88   2.06  1.94 2.63   2.37  2.61  1.81  1.31  1.21 1.89  1.67  1.95  1.09  2.17  0.861
 chr1   96493  96796 0.976 0.325 0.943 0.216  1.33  1.77 1.66   1.33  1.47  2.39  1.76  2.37 0.939 0.198 1.37  0.439 0.868 0.756
 chr1  271114 271432 1.95  2.62  2.04  2.29   2.49  2.86 2.12   2.89  2.30  2.25  3.16  3.50 0.695 1.03  0.491 0.785 1.15  1.00
 chr1  273026 273333 0.740 1.88  1.48  1.71   2.50  2.88 2.81   2.66  2.89  2.07  3.02  3.61 1.24  1.40  0.984 1.35  0.651 1.53
 chr1  274265 274599 0.954 1.59  2.06  0.227  2.05  1.19 0.906  1.60  1.62  1.99  2.05  2.47 0.920 0.211 1.06  0.426 1.99  0.529
 chr1  402478 402771 2.98  3.26  3.01  3.07   2.51  2.76 3.29   2.65  2.34  3.03  3.00  2.99 2.54  3.04  2.88  2.96  2.94  2.68

```

I tried inner join ,merge etc but i can map the chromosome cooridnates but the expression values comes empty

```
tmp <- merge(a,Rlog_ATAC_location, by=c("Chr","Start" ,"End"), all.x=TRUE)

Chr  Start    End CMP1 CMP2 CMP3 CMP4 CMP5 CMP6 GMP1 GMP2 GMP3 GMP4 GMP5 GMP6 HSC1 HSC2 HSC3 HSC4 HSC5 HSC6 Mono1 Mono2 Mono3
 chr1  16104  16348   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1  96494  96796   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 271115 271432   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 273027 273333   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 274266 274599   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA
 chr1 402479 402771   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA   NA    NA    NA    NA

```

Tried left join as well i get same result.

I want to map the coordinates as well as the samples columns which contains the expression

Any suggestion or help would be really appreciated

# Answer
> 2 votes
In your ATAC data frame, you have `chr1 16103 16348`. But in your merged data frame, you have `chr1 16104 16348`. Is it possible that you `a` is 1-based whereas your ATAC data frame is 0-based? This can happen in bioinformatics and can be annoying. This might be why you have a lot of NA in your merged data frame. You should subtract 1 from the start position of your `a` data frame while keeping the end position unchanged. Then try to merge again.

You should always be aware of your data format to know whether it is 1-based or 0-based.

---
Tags: r
---


title: PDB format: remark number for free text

text:
# Question
Title: PDB format: remark number for free text
I would like to add a text to PDB files that I'm processing with my tool, rna-pdb-tools. Someone points that the way I'm using it right now it's not correct (https://github.com/mmagnus/rna-pdb-tools/issues/48).

I use HEADER right now which is not correct.

```
HEADER Generated with rna-pdb-tools
HEADER ver 37c5b4e-dirty
HEADER https://github.com/mmagnus/rna-pdb-tools
HEADER Mon Oct 10 22:18:19 2016
ATOM      1  P     C A   1     -19.687  -3.296  65.469  1.00  0.00           P

```

Do you happen to know which remark number to use? (http://www.wwpdb.org/documentation/file-format-content/format33/remarks1.html)?

# Answer
From this document, it looks like REMARK 250 is the way to go.
> 4 votes

# Answer
**EDIT:**

Actually, after having a more thorough look into the documentation; `REMARK 3` might be the more relevant to what you need. Other tools are reported in that section.

> REMARK 3 presents information on refinement program(s) used and related statistics. For non-diffraction studies, REMARK 3 is used to describe any refinement done, but its format is mostly free text.

More on `Remark 3`

---

From the same source as you reference yourself v 3.3, I would argue that `REMARK 0` would be more accurate than `REMARK 250`.

> REMARK 0 (updated), Re-refinement notice REMARK 0 identifies entries in which a re-refinement has been performed using the data from an existing entry. This remark also describes the PDB code and the journal records for the original data set.

While when you consider the definition of `REMARK 250`

-- which follows directly:

* `REMARK 205` (specific to Fiber diffraction experiment),
* `REMARK 201`(and `REMARK 215/217` all specific to NMR experiment),
* `REMARK 230` (neutron diffraction study),
* `REMARK 240`(electron crystallography study) and
* `REMARK 245`(and `REMARK 247`: specific to EM study)

> REMARK 250, Other Type of Experiment Details
>
> REMARKs specific to other kinds of studies, not listed above. REMARK 250 is mandatory if other than X-ray, NMR, neutron, or electron study. The format of the date in this remark is DD-MMM-YY. DD is the day of the month (a number 01 through 31), MMM is the English 3-letter abbreviation for the month, and YY is the year.

and then `REMARK 265` is also about the crystallography experiment.

Hence, I know that personally, I usually assume that whatever is annotated in `REMARK 250`is relevant to the first acquisition experiment to me.

Another reason I am more inclined to think `REMARK 0`is more accurate, it is due that the definition of `REMARK 0` was **updated** while `REMARK 6-99` are *are no longer for use of free text annotation*.

BTW, I am not sure how Pdb\_extract would do with the remarks when converting the PDB file to CIF files. It would probably worth perusing PDB to PDBx/mmCIF Data Item Correspondences page.
> 2 votes

# Answer
**Use any REMARK number not reserved for something else.**

The REMARK field was initially for general remarks. In the format description from 1992 only the numbers 1 to 3 were reserved for specific things, with the note: *Other general commentary is given in higher numbered REMARKs*.

Then, in the next years, other remark numbers got reserved (and sometimes retired) for various specific annotations. The content of remarks was getting more and more formalized, to make it easier to parse by computers. The specifications from 1996 and 1998 recommended:

> Non-standard remark annotations, or those with no clearly-defined topic or assigned remark number, appear with remark number 6 or greater, but less than remark number 100.

I'd still follow this advice. It is generally understood that when parsing PDB files one may have custom REMARK records. I've seen also programs that output remarks without numbers - and this also doesn't cause problems with any reasonable PDB reader. On the other hand, using remarks with numbers reserved for very specific things (such as remark 250 in the most upvoted answer) can cause confusion.

The current spec (from 2012) is a precise description of files served by the wwPDB. It doesn't attempt to define a format for other uses, such as passing data between programs. In the section about REMARKs 6-99 it has the last bullet point saying: "REMARKs 6-99 are no longer for use of free text annotation", because these remarks were removed from the old PDB entries. But it doesn't prevent from using it outside of the wwPDB.
> 1 votes

---
Tags: file-formats, pdb
---


# Question
Title: Understanding some of the computational bottlenecks of Covid-19 research
I am a researcher in high-performance computing (with very little bioinformatics background), and I am trying to understand what are the current biggest computational bottlenecks of software used for research on the topic of Covid-19 (testing, vaccine development, drug identification, etc...).

I know this is a very broad question, but I am looking for the biggest common denominators. Said differently, if you could pick a few operations / algorithms / data structures to speed up Covid-19 research and they would magically run 1000 times faster, what would they be?

Note: links to books, papers, and scientific literature are very welcomed.

# Answer
I can only speak of drug design (and even then I am terrible at turning down the jargon).

In the case of drug design, this is pretty much plan C. Namely, none of compounds that entered clinical trial at the start of the year work (let's call this plan A) and none of the vaccines that are entering now clinical trial work (let's call this plan B although vaccination prevents so is a better solution). Clinical trials take a lot time as they need to be thorough (cf. the infamous thalidomide from the days with more lax testing)

There are several bottlenecks in drug design. These are technical based on our imperfect understanding of things, which results in less than ideal models. CPU time is not really an issue. Empirical validation is required at many steps along the way. This is greatly increased by automation. The XChem team at Diamond light source have done a fragment-screening against the SARS-CoV-2 main protease at an insane pace: things are changing fast!

## Drug screen

There are many ways to make a drug. In one way, a test system was set up (a reporter assay) that gives you a measurable read-out, say mm of clearance of bacteria on a Petri dish (viability) around a cellulose dish soaked in different separated natural extracts —and you might find that *Penicillium chrysogenum* has a compound of interest. What you are doing is drug screen. Screening a large panel of compounds (of varying purity) does not necessarily require knowing the protein target.

## Rational drug design

Fast forward to this decade and it becomes possible to predict the free energy release upon binding of a small molecule against a protein —mostly.

This requires a solved 3D structure of the protein by X-ray crystallography. Crystallising a protein is a tricky business, as a pure protein sample has to come out of solution in a folded form to form a neat lattice as opposed to just making unfolded gunk. What salt conditions work? Anyone's guess. Some have done deep learning on it, but the best solution is to have a robot try a very large range of combinations.

Once you have a protein structure you can predict what fits into the part that you want to target, say an active site or an interface. This is not done by simple cartesian geometric shape fitting like LEGO bricks as there are many molecular forces involved. So a mathematical model of these forces is made, called a force field. Many have been made, a lot of study has gone into these, but they are not perfect. It is all made worse by the fact that we cannot really model water properly —most models use a proxy (implicit solvent) or use non-polarisable water. This is all while made worse by the fact that a protein is not static but wobbles (following the same mathematical models). So a ligand does not slot in really like a key in a lock, but has an induced fit. Modelling all this is accurately is problematic, but one can still try: calculating where exactly a ligand binds to a protein is called for which one gets energy difference upon binding. And one can try an inordinate amount of small molecules —Zinc DB contains 230 million purchasable compounds. This is called a virtual screen. One can even start with the natural substrate and "weaponise it" by adding chemical groups that do a chemical reaction that breaks the protein. However, a metabolite is very different from a drug (cf. Lipinski's rule of five) as the latter must cross membranes etc. But generally the accuracy is not great without multiple steps of empirical validation.

To improve the results from a virtual screen the most common approach is to do a fragment screen to guide it. This is to automatedly soak protein crystals with a large library of small compounds and see which crystal contains a bound hit. The hits can be expanded into bigger fragments in a virtual screens and the short list of compounds tested (these are not cheap) by soaking and crystallising again. In the case of Covid19 virus, the cool thing is that this whole process is open in the aforementioned current campaign (Diamond is a government owned synchrotron and not a big pharma). What is really nice about this is that for the second cycle comp-chemists from everywhere have done their own guided virtual screen and have submit compounds (which you can see publicly).

It should be said that drug screens that are not rational design are still a thing. The recent antibiotic-discovering deep-learning paper in Cell\] is such an example, although it should be said that it was a feat not only in AI but also in automated empirical testing.

## Future

**Automation** via better and open software will definitely improve. With a small drive possibly from the dissemination of easy to use and to code off-the-shelf electronics, such as Rasperry Pis or Arduinos. But mainly with more complex machines such as the Opentron system (a python-controlled robot).

But indubitably where the greatest improvement can come is **better mathematical models**. Possibly even with deep learning. After all researches were amused last year-ish (2018 competition) when DeepMind "won" the protein 3D structure (fold) modelling competition (CASP13) and there have definitely been some really nice results in the field of docking with hybrid models. So deep-learning may still surprise us yet...

Lastly, **open science** will be a strong driver: using premade packages or tools on GitHub and code snippets from StackOverflow etc. will accelerate these.
> 5 votes

# Answer
**Drug screening**

I understand the idea behind this, particularly with deep learning playing such an important role in computation (see GPU). Moreover, molecular dynamics is famous for sumping computer power, so yeah this is a bottleneck. However, this is not trivial because the wet-lab work is such a long way off. The first crystal/cyro/EM/NMR structure that will be resolved is **Spike** (S) and this would not be a 'traditional' drug target, its a vaccine candidate (okay influenza has surface antigen drug target, Tamiflu). Expressing S is a priority (vaccines) and solving the structure is a by-product and its not that big, plus the majority of it is water-soluable so it will cystalise nicely (its a surface antigen). Obtaining the structure for a more 'traditional drug target' which would be ORFa and ORFb (RNA dependent RNA polymerase ... and other replication associated proteins), will be difficult, because they are huge (ORFa + ORFb are 20 Kb of RNA, which for an RNA virus is twice as big as the genomes for most RNA viruses) and for crystals ... solubility(?). Ultimately, how do you cystalise a protein that 4kb of amino acids - well you'd have to do it in subdomains, you'd probably need NMR/ really hot EM.

**General**

HPC is a big field and it highlights the importance of one area most struggle with which is GPU (graphic processing unit) for deep learning (DL) applications. HPC's role molecular dynamics/drug discovery) and a mixture of DL and MD.

Traditional CPU 'farms' most of use have access to, is less critical because whether it takes 2 days per calculation or 2 hours isn't a huge issue (but its very seriously annoying indeed)..

The most pressing need for computational power will be epidemiological modelling, particularly if someone is pursing a microsimulation. The firepower would enable them to explore a large number of scenarios and parameter explorations. Apart from a vaccines then epi modelling is the hot-topic of the moment.
> 1 votes

---
Tags: sequence-alignment, sequence-analysis, software-recommendation, covid-19, sars-cov-2
---


# Question
Title: Are there any databases of templates for common bioinformatic file formats?
I want some templates of different file formats that I can use to test my scripts and identify possible bugs in my code.

For example, consider nucleotide FASTA, a simple but often abused format, I would want templates to capture regular and irregular formats, like I have seen all of these:

1) Single line sequence

```
>1
ATG

```

2) Multi-line sequence

```
>1
AT
G

```

3) Upper and lower case letters in sequence

```
>1
Atg

```

4) Ns and Xs (and possibly other letters) in sequence

```
>1
ANnxX

```

5) Unusual headers (sometimes non-ASCI characters, need to consider the encoding)

```
>ATG >汉字
ATG

```

6) Whitespace between records

```
>1
ATG

>2
ATG

```

7) Duplicated headers

```
>1
ATG
>1
ATC

```

8) Empty headers or sequences (valid FASTA?)

```
>
>

```

9) No new line `'\n'` character on last line (can mess-up file concatenation)

```
>1
A# < no new-line here

```

10) Different newline characters depending on the OS

```
>1
A# \r\n vs \n

```

etc.

There should be separate templates for nucleotide and protein FASTA, and separate ones for aligned FASTA.

It would ideally include other aspects too, like different compression formats (such as `.gz`, `.bzip2`) and different file extensions (such as `.fa`, `.fasta`).

I have never seen resources that provide templates covering these, but I think it would be useful. Of course I could build my own templates but it would take time to capture all the likely variations of the formats, particularly for more complex file formats.

Note, I am not just interested in FASTA format, it was an example.

Also note, I know about tools (such as `BioPython`) that should handle many formats well, but they may also have bugs. Anyway, in practice sometimes I end up parsing files myself directly because I don't want the overhead or dependency of an external package.

*EDIT: Please don't answer this question to say you don't know of any such resources, me neither, hence the question. bli's helpful answer shows that there is at least one test suite that could be used as a starting point. I know that it is normally easy to look up the specification of any particular file format.*

# Answer
> 2 votes
BioJulia maintains BioFmtSpecimens. From the README:

> Bioinformatics is rife with formats and parsers for those formats. These parsers don't always agree on the definitions of these formats, since many lack any sort of formal standard.
>
> This repository aims to consolidate a collection of format specimens to create a unified set of data with which to test software against. Testing against the same cases is a first step towards agreeing on the details and edge cases of a format.

# Answer
> 7 votes
You mention Biopython, which contains tests: https://github.com/biopython/biopython/tree/master/Tests.

Some of the tests consist in reading files present in the folders listed in the above link. These files could be a starting point for a database of test files. Whenever one comes across a test case not covered with these files, one could construct a new test file and contribute it to Biopython, along with a test, or at least file an issue: https://github.com/biopython/biopython/issues.

That would be a way to contribute to Biopython while constituting a database of test files.

# Answer
> 5 votes
Not that I am aware. It is best to go with format specifications when coding.

Also it may be good to look at the example files that come together with various tools performing file conversions and handling. E.g.

# Answer
> 2 votes
As far as I know, there is no single repository that collects all of the common data formats used in bioinformatics. Typically, you have to go to the source to find the specifications for each format. There are a few places that collect descriptions of file formats, though:

* IGV File Formats, coveres all of the formats usable in the Broad Institute's Integrative Genomics Viewer software (which is a lot).
* NCI File Formats, mostly formats used by TGCA (including MAF and VCF).
* UCSC Genomics, covers BED, MAF, and a few others.
* GenePattern, covers many of the file formats related to microarray data.
* GSEA, Broad Gene Set Enrichment Analysis documentation.

# Answer
> 1 votes
No. At least none that I've heard of and I doubt there ever will be. There is no central repository for formats and each tool, community, field etc have their own.

The best you can do is look up the official standard for each format and hope they include examples. Having a truly comprehensive collection of all possible variations is basically impossible for ~~horrible~~ complex formats like VCF. Just consider the simple fact that it allows user-defined `INFO` fields with arbitrary contents.

So the best you can do is make sure your scripts conform with the standard and hope your input does as well.

# Answer
> 0 votes
These may not be exactly what your are looking for, but they do contain a wide range of formats with examples. If you want more, you might try searching other major databases hosting other types of into e.g. Uniprot, PDB, NCBI.

https://genome.ucsc.edu/FAQ/FAQformat.html

http://www.ensembl.org/info/website/upload/bed.html

---
Tags: fasta, database, file-formats
---


# Question
Title: Cuffmerge: EOF marker is absent. The input is probably truncated
I need to merge my all transcripts.gtf from cufflinks output follow this command line :

```
cuffmerge -o merged_gtf_output -p 15 -s ref.fasta -g anot.gtf assembly.txt

```

But I have this warning :

```
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
File merge_B3//tmp/mergeSam_fileOaNffV doesn't appear to be a valid BAM file, trying SAM...
[21:40:04] Loading reference annotation.
[21:40:05] Inspecting reads and determining fragment length distribution.
Processed 3921 loci.
> Map Properties:
>       Normalized Map Mass: 33242.00
>       Raw Map Mass: 33242.00
>       Fragment Length Distribution: Truncated Gaussian (default)
>                     Default Mean: 200
>                  Default Std Dev: 80

```

Can someone tell me what this warning means and how to fix it?

# Answer
> 1 votes
There is a post about this on seqanswers. The poster says the following about cuffmerge behavior:

> I kept the temporary files while running cuffmerge. It builgs a SAM file from transcripts.gtf files. Then it seems first it reads it as BAM since I read samtools take the files as BAM default. Then it finds it is not a BAM so it tries SAM. Then it can read it and resume the task.

So it sounds like the BAM file is generated from GTF as a kludge inside cuffmerge, not the other way around.

Did you change any of the intermediate files or delete anything? Is there maybe a different cuffmerge run that it might be colliding with?

Did you get any output? it sounds like this error is frequently observed when actually the command completes fine, e.g. https://www.biostars.org/p/189139/. See other posts by googling the problem maybe, e.g. https://www.biostars.org/p/344439/

(I know this isn't a proper answer but it was too involved for a comment.)

---
Tags: rna-seq, samtools
---


# Question
Title: Is there a aminoacid mutation that forms similar protein?
Is there an amino-acid shift that will 99% of the time end up resulting the same biochemical function/structure and implications? Like Alanine replaced by Leucine will code in 99% of the cases(at any position that first has Alanine) towards proteins that will function the same without big percentual changes in how good it binds or catalyses.

I do have a big amount of mutations. I'd like to filter them out on chance of biological/phenotypical impact, and since some amino acids look like good substitutes by eye (if I look at the molecular structure)

# Answer
A mutation is most often likely to be either neutral (no effect) or destabiling (protein misfolds). Changes in activity require changes at the active site or the entranceway to it —which are often only a handful of residues, so are very uncommon.

So the answer to your question is most mutations keep the activity, although a large fraction break it. While to get a new activity the process is more fluid and is not boolean.

NB. Not all proteins are enzymes, but for simplicity I will concentrate on enzymes.

## Neutrality and core vs. surface

The outer residues of a protein, the surface residues, are very tolerant to mutations. Unless they are part of an binding surface (interface) or the change results in a charge change that is. The core residues on the other hand are more sensitive to mutations.

## Destabilisation

One can do a scan *in silico* or *in vitro* or *in vivo* to determine the mutational landscape of a protein, which can take various forms.

This being Bioinformatics stack overflow I'll give an in silico example. The following for example is a heatmap I did some time back with Rosetta pmut\_scan, which is a classic tool for this and shows in the y axis the residue it was mutated to and on the x axis the position. A larger change in Gibbs free energy of folding) relative to wild type (∆∆G) means it is destabiling, but anything less than 5 kcal/mol is neutral. As you can see many are neutral, but some half are destabiling —with a caveat for proline, which require a backbone change.

In the in vitro setting one speaks of melting temperature (DSC assays) or turnover number over Michaelis constant (enzyme specific assay). In the in vivo setting one speaks of a fitness cost of a mutation, i.e. the protein is less active so the organism grows less.

## Drug resistance

Drugs inhibit enzymes by binding the active site. Therefore resistance is often costly. There are some nice papers that look into the fitness cost and mutational landscape of HIV protease or β-lactamase. In this situation a strong epistasis is seen. That is, a first mutation that confers resistance to the drug greatly destabilises the protein, while a second compensates for the first mutation.

## Change in function

An enzyme can catalyse a reaction well on a given substrate, because it was evolved this way (technically the encoding gene evolved, but shh). This is the physiological activity. However, enzymes can also catalyse side activities on similar substrates with a lower turnover and higher Michaelis constant. This is especially true for substrates it was not selected against. This is called "enzyme promiscuity". The tolerance of the physiological activity to mutation is called "robustness", while the ability to increase promiscuous activities is "plasticity". Different enzymes have a different balance of the two.

So a mutation may change the strength of the physiological activity, but may result in a new promiscuous activity that can be used as a stepping stone for selection to increase this activity. Therefore, as you can see it is not a binary switch, except for a few very specific case.

## Best amino acid for changes?

Proline is a very strange amino acid as it is a secondary amino acid. It lacks a hydrogen on the backbone nitrogen so it cannot form the required hydrogen bond for an α-helix. And can adopt a cis-conformation with minimal penalty (whereas it is about 5 kcal/mol for other chiral amino acids). So changes to this amino acid are often very destabilising.

Protein cores tend to be hydrophobic so changes to the hydrophobic residues are often deleterious, especially when going from a small residue to a larger one.

The remaining residues are frequently found on the surface so may be neutral. However, some are found in the active site and may be able to shift the parameters of the ligand space for that enzyme. A change in charge may result in an increased specificity for a different set of substrates for example. A size change from a larger amino acid to a smaller one may result in a broader range of promiscuous activities —mutations to alanine for example.

So ironically, polar and charged residues as they are frequently found on surfaces are the residues which are most often neutral, but in a small amount of cases they may change activity.
> 4 votes

---
Tags: mutations
---


# Question
Title: nextflow `fromPath` not passing all files with *
main.nf is

```
#!/usr/bin/env nextflow

params.input = ''
tchan = Channel.fromFile(params.input)

tchan.println()

```

I have a.bam, b.bam and c.bam

I run:

```
nextflow run main.nf --input *bam

```

But it only prints `a.bam`

According to the documentation, wildcard character it is interpreted as a glob path matcher

What am I doing wrong?

Thank you

# Answer
got it!

I should run:

```
nextflow run main.nf --input "*bam"

```

instead of

```
nextflow run main.nf --input *bam

```
> 2 votes

---
Tags: nextflow
---


# Question
Title: truncate fastq.gz file without error in Snakemake
All,

As the first step in my Snakemake pipeline, I generate small fastq files from all input fastq.gz files. This is to run a quick round of "infer\_experiment.py" (from RSeQC, to get the kit strandedness). However, Snakemake says this step fails. I think I know why, it is because a command like this

```
zcat fastq_raw/20190011-043_S43_R1_001.fastq.gz | head -n 400000 > infer_strandedness/20190011-043_S43_subset_R1.fastq

```

Ends with an error: `gzip: stdout: Broken pipe` and Snakemake uses Bash strict mode. (Snakemake reports: `(one of the commands exited with non-zero exit code; note that snakemake uses bash strict mode!)`). So my question is:

How to avoid this, or, are there other truncating methods that do finish cleanly? Preferebly they are also efficient (seqtk offers subsampling but this parses the entire fastq file!).

# Answer
Snakemake is using bash strict mode to ensure best practice error reporting in shell commands. Unfortunately some tools violate bash strict mode.

More details of the error, here. Alternatively, you can use if you don't want disable the strict mode.

```
zcat file.fastq.gz  | awk '(NR<=400000)' > subset.fastq

```
> 1 votes

---
Tags: snakemake
---


# Question
Title: How did researchers derive the Ramachandran "validation" contours?
I'm a beginner to structural biology and for fun, calculated the tortional angles of some 100, 000 proteins. Here is my Ramachandran plot:

When I went to look for "canonical" Ramachandran plots, I discovered most researchers overlay a "contour" plot of allowed torsional angles (see below). My two questions are:

* How are the distributions (ie. **countours**) of protein secondary structures calculated/derived? (ie. How do I "cluster" these into beta-sheets, alpha helices?)
* How can I derive them myself?

# Answer
To expand on Matteo Ferla's comment as an answer, there are two approaches to getting these "allowed" contours:

1. From probabilities found in observed crystal structures. You could imagine plotting a Ramachandran plot for the whole PDB and colouring the points by assigned secondary structure from DSSP. The contours would then be regions that encompass e.g. 99% of alpha-helical residues, or 99% of all residues. In doing this you would notice for example that glycine can adopt more regions of the Ramachandran space that other amino acids. PROCHECK is software that scores a protein based partly on the likelihood of the phi/psi angles of each residue.
2. From reasoning about the chemistry of the secondary structure of proteins. Only certain phi/psi combinations allow the i -\> i + 4 pattern in an alpha-helix or the alternating up/down sidechain pattern in a beta-strand.
> 3 votes

---
Tags: structural-biology
---


title: Running a script in python

text:
# Question
Title: Running a script in python
I am trying to run a script in python but I getting this error

```
loan-mac-13:vcf2maf fi1d18$ python
Python 2.7.10 (default, Feb 22 2019, 21:55:15)
[GCC 4.2.1 Compatible Apple LLVM 10.0.1 (clang-1001.0.37.14)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import os
>>> os.chdir("/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf")
>>>
>>>
>>> os.getcwd()
'/Users/fi1d18/Downloads/vcf2maf-master/vcf2maf'
>>> vcf2maf.py -h
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
NameError: name 'vcf2maf' is not defined
>>> python vcf2maf.py -h
  File "<stdin>", line 1
    python vcf2maf.py -h
                 ^
SyntaxError: invalid syntax
>>>

```

Where I am doing wrong?

Thanks for any help

By your comments I run

Sorry says that Last login: Wed Apr 1 11:50:14 on ttys000

```
pc-133-235:~ fi1d18$ pwd
/Users/fi1d18
pc-133-235:~ fi1d18$
pc-133-235:~ fi1d18$ cd Downloads
pc-133-235:Downloads fi1d18$ python vcf2maf.py -h
Traceback (most recent call last):
  File "vcf2maf.py", line 19, in <module>
    import vcf
ImportError: No module named vcf
pc-133-235:Downloads fi1d18$

```

# Answer
You are trying to run a script file in the python console. That doesn't make sense. The console is for running commands manually, interactively. A script is a collection of commands. You can't just run that by name in the console. If you want to use it, copy the *contents* of your script and try to run them in the console.

The way you run a python *script file* is simply:

```
python vcf2maf.py -h

```

Or, assuming you have made the file executable and it has the right shebang line, just:

```
./vcf2maf.py -h

```

---

The specific error you are getting is because `foo.bar` in Python means "run the `bar` method from the object `foo`. So `vcf2maf.py` means run `py` from `vcf2maf`. Since there is no object defined with the name `vcf2maf`, you get the `name 'vcf2maf' is not defined`.
> 4 votes

---
Tags: python, vcf, software-usage
---


title: Allele count is zero

text:
# Question
Title: Allele count is zero
I ran freebayes on bacterial data and got a vcf file. Some of the SNPs found have `AC=0`, but in the bam file I can see that there are alternative alleles, so in my understanding it supposed to be `AC>0`.

example: in vcf file:

```
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  55509   55508   55506   55507
scaffold1       78      .       A       C       9.27683e-13     .       AB=0;ABP=0;AC=0;AF=0;AN=3;AO=466;CIGAR=1X;DP=2254;DPB=2254;DPRA=0;EPP=20.9226;EPPR=29.3274;GTI=0;LEN=1;MEANALT=2;MQM=50.3412;MQMR=45.5379;NS=3;NUMALT=1;ODDS=2080.75;PAIRED=0.899142;PAIREDR=0.533932;PAO=0;PQA=0;PQR=0;PRO=0;QA=16347;QR=63107;RO=1783;RPL=203;RPP=19.7856;RPPR=291.853;RPR=263;RUN=1;SAF=465;SAP=1006.25;SAR=1;SRF=894;SRP=3.04075;SRR=889;TYPE=snp       GT:GQ:DP:AD:RO:QR:AO:QA:GL      .       0:126.704:776:626,149:626:22118:149:5255:0,-1217.7      0:126.704:834:663,170:663:23671:170:6004:0,-1224.79     0:126.704:644:494,147:494:17318:147:5088:0,-770.853

```

but with samtools tview:

Can someone tell me if this makes sense and why it happens?

# Answer
> 3 votes
The AC is the number of alternate alleles from called genotypes, you can see this in the vcf header:

```
##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of
alternate alleles in called genotypes

```

Assuming diploid, if we have 3 individuals, 0/0 , 0/1, 0/0, AC will be 1, because it detects 1 count of alternate allele.

If you look at your genotype calls, all are 0, so AC = 0

---
Tags: vcf, snp
---


# Question
Title: protein secondary structure matching/comparison scoring
I want to compare two secondary structures of aligned proteins. I do not know exactly how to do it well.

Example:

```
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC

EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHEEEECCCCCCCCCCEEEEEEEEEEEEEE

```

Example of matching:

```
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
|||        |||||| |||                    |||||||||     ||||||      |   |||    ||
EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH----HHEEEECCCCCCCCCCEEEEEEEEEE

```

Hypothetically: 50% aligned - \> score 0.5 - \> thats scoring is wrong:

```
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
||||||||||||||||||||||||||||||||||||||||||
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHEEEEEEEEEEEEEEEHHHHHHHHHHHHHHHHHHHHHHHH

```

This case will have similar score...

This is a puzzle for me, I don't know how to approach it technically.

# Answer
Hey i solve this simply by use of SOV alghoritm. https://www.researchgate.net/publication/319526145\_A\_modified\_definition\_of\_Sov\_a\_segment-based\_measure\_for\_protein\_secondary\_structure\_prediction\_assessment

cheers
> 1 votes

---
Tags: proteins, python, statistics, secondary-structure
---


title: How can I compare two bed files?

text:
# Question
Title: How can I compare two bed files?
I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's LiftOver tool and CrossMap but saw that they give me different results. I therefore need a way of assessing how correct the results are.

I decided to test by downloading a bed file of the RefSeq genes for hg37 and one for hg38 from UCSC. Then, I run `liftOver` to translate the hg37 file to hg39 coordinates and now I want to compare the results of the liftover to the bed file I downloaded for hg38.

How can I usefully compare these files? I am hoping for a tool that can report how similar the two are. Ideally, I would like to see three things in the output:

1. Matching regions: those that are shown with identical positions in the liftover results and the UCSC bed file.
2. Non-matching regions: those whose coordinates do not match between the two files.
3. (if possible) fuzzy matches: those that are not identical but pretty close (off by a few nucleotides).

I thought it would be easy to do with a little awk script but it gets quite complicated because certain RefSeq genes are mapped to multiple positions (e.g. MIR4454). And the 3rd requirement (fuzzy) is not trivial to do in awk. I can live without that though, if a tool can give me the other two.

# Answer
> 3 votes
From a set operations viewpoint, consider your hg37 regions as a "reference" set and hg38 as a "map" set.

If you have these sets in BED format, BEDOPS lets you do set operations on them with various overlap criteria.

Specifically, you can use BEDOPS `bedmap` to map hg38 regions to hg37 regions. The `bedmap` tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap.

Once you do mapping, you can do operations, like counting the number of regions that meet the overlap criteria between reference and map sets, and calculate percentages.

For example, use `--exact` to get exactly-matching regions between hg37 and hg38:

```
$ bedmap --count --exact hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt

```

You can do fractional ("fuzzy") overlaps, where (for example) at least 50% of an hg38 region has to overlap an hg37 region:

```
$ bedmap --count --fraction-map 0.5 hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > fuzzy_0p5_counts.txt

```

Getting a percentage out of a "fuzzy" result can be complicated by overlaps between hg38 regions that overlap a reference hg37 region. A more complicated procedure can be implemented with the `--echo-map` option and some command-line work to implement logic to deal with that case more stringently.

To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using Kent tools `fetchChromSizes` and BEDOPS `bedops --difference`:

```
$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - > hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed > hg38.masked.bed

```

The regions in `hg38.masked.bed` make up all of the genomic space in hg38, except for your lift-over regions.

Apply `--count --fraction-ref 1.0` to count the number of reference (hg37) elements which completely overlap these masked hg38 regions:

```
$ bedmap --count --fraction-ref 1.0 hg37.bed hg38.masked.bed | awk '{s+=$1}END{print s;}' > no_counts.txt

```

The documentation for `bedmap` describes these overlap criteria and operations in more detail.

# Answer
> 2 votes
I wrote ParsEval to handle these types of comparisons. ParsEval reports a variety of similarity statistics at both the nucleotide level and feature (exon) level. It doesn't explicitly support "fuzzy" matches, but it does report scaled values between 0.0 and 1.0 for similarity, so you could get a similar feel by definig a similarity cutoff for filtering.

Caveat: ParsEval doesn't handle BED files, only GFF3.

# Answer
> 0 votes
I was trying the code from @Alex Reynolds, however, I got retrieved an error:

```
bedmap --count --exact SRR800856_*.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt
May use bedmap --help for more help.

Error: Option SRR800856_ngs_te_mapper_nonredundant.bed does not start with '--'

```

But when I tried with just two bed files:

```
bedmap --count --exact SRR800856_ngs_te_mapper_nonredundant.bed SRR800856_popoolationte_nonredundant.bed

```

I don't get any error. Could you tell me what I should add. Thank in advance for your time. P.D. As I am new here I couldn't get to comment on the specific thread.

---
Tags: bed, awk, liftover
---


title: qvalue and p.adjust functions

text:
# Question
Title: qvalue and p.adjust functions
Between the qvalue package and the p.adjust function, which is more appropriate to use when trying to calculate the q-values of a dataset? According to the manual for the q-value package, the q-value calculated is not an adjusted p-value which is what the p.adjust function would return. There seems to be an equal amount of people that use one or the other, but I wanted to see other opinions that could help clarify what would be most appropriate to use.

# Answer
> 0 votes
`p.adjust` can calculate adjusted p-values with different methods. I think you can choose Bonferroni correction and BH correction, for instance. BH correction gives you False Discovery Rate (FDR), which is one of the most commonly used I guess.

Q-value is different and cannot be calculated by `p.adjust`. The idea is that when you have many p-values to adjust, they form a distribution. If all p-values come from the null hypothesis, the p-value distribution is going to be uniform on \[0, 1\]. But usually, our p-values are a mixture of the null distribution and alternative distribution. The parameter $\pi\_0$ describes the proportion of p-values (or tests) that comes from the null hypothesis. To calculated q-values, the algorithm first estimates $\pi\_0$, then calculates the q-value.

Note that if you set $\pi\_0=1$, which means that you believe all the p-values come from the null distribution, then q-value is equal to FDR.

I think most people do not know about the meaning and details of q-values and it might be sufficient to use BH procedure in most cases.

---
Tags: r, bioconductor
---


# Question
Title: How can I obtain the percentage gene expression per identity class in Seurat as further processible numbers (e.g. matrix)?
I am analysing my single cell RNA seq data with the Seurat package. I want to know if there is a possibilty to obtain the percentage expression of a list of genes per identity class, as actual numbers (e.g. a matrix) which I can write out to say an excel file. For example I want to be able to say gene X is expressed in 30% of cell A, 10% of cell B etc.

I can generate something comparable with:

```
DotPlot(object = my_object, genes.plot = "my_gene")

```

However the results are only graphic and I wish to have further processible numbers.

Furthermore:

```
AverageExpression(object, genes.use = "my_gene")

```

Produces expression values which I cannot transform to percentages.

I will be very grateful on any hints.

# Answer
This can be solved like this:

```
library(Seurat)

my_genes <- c("gene1", "gene2", "gene3")

exp <- FetchData(object, my_genes)

matrix <- as.matrix(colMeans(exp  > 0))*100

```
> 5 votes

# Answer
You can get the table that is used to make the dot plot if you modify the DotPlot function to return it instead of the ggplot, and use the argument do.return=T.

To edit the function, the command is:

```
trace("DotPlot",edit=TRUE)

```

Then replace the last line "return(p)" by "return(data.to.plot)" and save the edit.

You can call the function and store the return in a variable.

The editing doesn't affect the source code, so the function will be untouched at your next R session.

Best,
> 3 votes

# Answer
Also potentially useful is the solution proposed here: https://github.com/satijalab/seurat/issues/1888

Basically, save DotPlot to an object, and your percentages are found within.
> 1 votes

---
Tags: r, scrnaseq, seurat
---


# Question
Title: What is the difference between fixed effects and random effects in the context of Linear-mixed models?
The terms in LME i.e. fixed and random create confusion? What is the genesis that can distinguish between the two ?

# Answer
I do not know much about statistics but I will try my best to explain.

First, random effects are defined as the factors (categories) in the population that we are not aware of (not observed), so we are randomly sampling levels of those factors when we sample the population.

Practically speaking, random effects can be found when there are hierarchical structures in the data. For example, we want to compare the relationship between students' SAT score and their college admission rates. You might just want to run a simple linear regression. But there might be some structure in the data that is hidden. For example, the students can be grouped into classes, and then schools. You can also group students by the social-economic status of their parents.

The idea is that students belonging to different categories may have different college admission rates even if they have the same SAT score. Therefore, it is important to consider some of these "random effects" when building your model.
> 1 votes

# Answer
Linear regression hypothesis is that effect of a variable is comprised as a constant component plus variability. In other words, the independent variable (s) is deemed to reflect a fixed (constant) part - determined as population effect size - and a variable part - random effect. Further, the term - Linear Mixed(regression)model is used to reflect use of continuous variables and nominal variables acting as independent variables that influence the dependent variable. THE CONTINUOUS VARIABLE is comprised of population effect size - and a variable part - random effect . The nominal or categories based variable effect is probabilistic and we use the term probabilistic variable for such variables. Sometimes academics mistakenly call its effect on dependent variable as random effect. Finally, this may be noted that principles of ANOVA and probability have been combined to formulate linear Mixed effects models.
> 0 votes

---
Tags: statistics, linear-regression, structural-biology
---


# Question
Title: Calculating the number of probes for a given genomic range
I have a big data frame like this for copy number (exome seq)

```
    > head(df)
  Chromosome    Start      End
1 1             64613  5707515
2 1           5712940  5732322
3 1           5732399 16383682
4 1          16383742 16389288
5 1          16390813 16830026
6 1          16830201 17278112

```

I also have coordinate of human exams like this

```
Exone   Choromosome     Start_of_exone  End_of_exone
uc001aaa.3  chr1    11873   12227
uc001aaa.3  chr1    12612   12721
uc001aaa.3  chr1    13220   14409
uc010nxr.1  chr1    11873   12227
uc010nxr.1  chr1    12645   12697
uc010nxr.1  chr1    13220   14409
uc010nxq.1  chr1    11873   12227
uc010nxq.1  chr1    12594   12721
uc010nxq.1  chr1    13402   14409

```

Actually I need to interset my copy number coordinates and human exome coordinates to see there are how many exomes in these coordinates

I have used this code but all returns 0

```
snp_table <- expanse coordinate file
genomic_ranges <- df

how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(snp_table$Chromosome == coords[1] & snp_table$Position > as.numeric(coords[2]) & snp_table$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_snps <- apply(genomic_ranges, 1, how_many_in_range)

```

Can you help me?

# Answer
> 0 votes
This command gave me the number of axons I hope this is correct

Please correct me if I am wrong

```
how_many_in_range <- function(coords){
    coords = as.vector(coords)
    sum(human_exome_coordinates$Choromosome == coords[1] & human_exome_coordinates$Position > as.numeric(coords[2]) & human_exome_coordinates$Position < as.numeric(coords[3]))
}

genomic_ranges$number_of_exon <- apply(segment_mean_table, 1, how_many_in_range)

```

# Answer
> 1 votes
Regardless of what your specific data is, the goal of your analysis is really to find overlaps between two data frames of genomic coordinations. My favorite package for doing this is GenomicRanges. It is very fast and can take care of different strands (if your data contains strand information).

---
Tags: r, bash, cnv, exome, genomicranges
---


title: A box plot of qualitative variables

text:
# Question
Title: A box plot of qualitative variables
I have 2 groups of patients; responders to chemotherapy and non-responders

I have calculated cancer cell fraction (CCF) of a set of driver genes (by variant allele frequency) for each group and I have

```
> head(dat)
    Response       CCF
1 Responders 1.0000000
2 Responders 0.5413323
3 Responders 1.0000000
4 Responders 1.0000000
5 Responders 1.0000000
6 Responders 1.0000000

unique(dat$Response)

[1] Responders    Nonresponders
Levels: Nonresponders Responders
>

```

If CCF \> 0.95 , the mutation is clonal otherwise sub clonal

I want to show how many clone and subclone are in these two groups by box plot or something similar like below

I have tried this which was nonsense

```
ggplot(dat, aes(x=Response, y=CCF)) +
    geom_boxplot()

```

Can you help me?

# Answer
I'm not sure that boxplot will be the more appropriate representation, as you will end with two numbers (count of Clonal and SubClonal) per groups of patients.

One solution will be to first create a new categorical variable based on CCF values using for example `ifelse` statement for example by writing:

```
df$Clonal <- ifelse(df$CCF > 0.95,"Clone","SubClonal")

```

Then you can get the count of each Clone and subclone for each responders and non responders by using `table`:

```
DF <- as.data.frame(table(df[,c("Response","Clonal")]))

```

You can finally convert these new values as percentage by grouping according to the Response. For example, using `dplyr`, you can do something like that:

```
library(dplyr)

DF %>% group_by(Response) %>% mutate(Freq_percent = Freq / sum(Freq))

```

And finally, you can plot it as a barchart or single point in `ggplot2`. Hope it helps you to figure it out how to deal with your data.
> 1 votes

---
Tags: r, ggplot2, software-usage
---


title: co-localisation error code using co-loc

text:
# Question
Title: co-localisation error code using co-loc
Can anyone offer advice on what the below error code relates to. I am using co-loc R package with 2 GWAS datasets. Using coloc.abf function with pvalues, MAF, sample size and case/control cohort. Thank you

```
  results <- coloc.abf(  dataset1=list(snp=as.character(1_subset[,2]),pvalues=1_subset[,11],N=1_data_subset$sample_size, type="cc",s=0.327,MAF=1_data_subset[,12]),   dataset2=list(snp=as.character(2_subset[,2]),pvalues=2_subset[,8],N=2_subset$sample_size,s=0.4319,MAF=1_subset[,12], type="cc"),p12=1e-06)

error code: Error in data.frame(V, z, r, lABF) :
  arguments imply differing number of rows: 1994, 1992
In addition: Warning messages:
1: In 2 * N * f :
  longer object length is not a multiple of shorter object length
2: In 2 * N * f * (1 - f) :
  longer object length is not a multiple of shorter object length
3: In r * z^2 :
  longer object length is not a multiple of shorter object length

```

# Answer
> 0 votes
In your `dataset2 = list()`, the MAF argument should be `MAF=2_subset[,12]` instead of `MAF=1_subset[,12]`. Since your two datasets have different numbers of lines, this is what causes the error.

I am not sure whether this is a typo in your code, or it is because you do not have MAF information in the second dataset. What I like to do before COLOC is to filter two datasets so they have the same number of lines, and then feed those filtered datasets to the coloc function.

(Also just out of curiosity, is `1_subset` a data frame for your GWAS results? I think it is not recommended to have variable names starting with a number.)

---
Tags: r
---


# Question
Title: merging two/or more bed file into one bed file
I am trying to merge two bed files (more in future) to one. my bed files are something like : .

I need to merge them in a way to have the shared chromosome location.

Is there a way to do that ?

# Answer
In my opinion, the easiest way to merge bed files is to use `bedtools merge`.
> 2 votes

# Answer
1. Convert each of your Excel spreadsheets to tab-delimited text files (how-to).
2. Install and run Cygwin, if using Windows. If you are using OS X, open up the Terminal app.
3. Install the BEDOPS toolkit to get sorting and set operation tools (if you're not averse to using them).
4. Strip out the Microsoft line endings, strip the header, and sort each of the *N* files:

```
$ tr -d '\r' < fileN.tsv | tail -n+2 - | sort-bed - > fileN.bed

```
5. Merge the files:

```
$ bedops --merge file1.bed file2.bed ... fileN.bed > answer.bed

```

The file `answer.bed` will contain all merged regions from files `file1.bed` through `fileN.bed`. Merging will be done on intervals for each unique chromosome name.
> 1 votes

---
Tags: sequence-alignment, samtools, sequence-analysis, bedtools, macs2
---


# Question
Title: run a tool with python>=3.6 on galaxy-project instance that runs on python2.7
I am developing a wrapper for a bioinformatics tool writen in python3.

The native Galaxy-project runs on python2.7.

How can I the tool that requires python3 in a galaxy-project instance?

# Answer
Each Galaxy tool uses its own (typically conda-based) environment, which will normally not use the same python as Galaxy itself (note that Galaxy can run under python 3). An example of this is deepTools, which is python 3 only but will happily run under a python 2 Galaxy since conda is being used to create a python 3 environment for it. It's best to distribute your tool via bioconda and then Galaxy will be easily able to install it.
> 2 votes

---
Tags: python, galaxy
---


# Question
Title: plot correlation between several bam files
I have several mab files from medip-seq. they supposed to be replicates and I have to draw their correlation. First I used multiBamSummary from deeptools to merge all the bam files based on bi/bed. Now when I am using the plotCorrelation I do not get any file. the argument is :

```
plotCorrelation --corData myMultiBamSumary.npz --corMethod spearman --whatToPlot scatterplot --skipZeros --removeOutliers --plotFile Myplot --plotFileFormat pdf

```

Is that the correct strategy to pot the correlation? Is there any suggestion?

Thanks alot

# Answer
> 1 votes
That's a completely normal way to plot correlation between them. Spearman correlation was a wise choice on your part, since it's rank based and you don't have to care as much about normalization.

Update: It turns out that matplotlib is broken and I need to have a new deepTools release to work around that. That should be available sometime Sunday.

---
Tags: sam, sequence-alignment, sequence-analysis, chip-seq, deeptools
---


# Question
Title: Getting the number of weeks from a time span
I have a very big clinical data like

```
Date.Death  Date Of Operation
02/08/2015  27/11/2012
02/08/2015  27/11/2012

```

I want to know how many weeks each patient has lived after operation

I then tried

```
> t1 <- as.POSIXct("27/11/2012")
> t2 <- as.POSIXct("02/08/2015")
> as.double(difftime(t2,t1,unit="weeks"))
[1] -1317.571

```

Or

```
> span <- interval(as.POSIXct("27/11/2012"), as.POSIXct("02/08/2015"))
> span
[1] 0027-11-20 LMT--0002-08-20 LMT
> t <- as.period(span, unit="day")
> t
[1] "-9223d 0H 0M 0S"
> t / as.period(dweeks(1))
estimate only: convert to intervals for accuracy
[1] -1317.571

```

Both give wrong number of weeks

Manually calculating one by one is painful and error prone

You know anyway to get the number of weeks?

Thanks

# Answer
I think your error is located in your use of `as.POSIXct` which gives you wrong date formats because you did not provide a format in your function:

```
> as.POSIXct("27/11/2012")
[1] "27-11-20 LMT"

> as.POSIXct("27/11/2012", format = "%d/%m/%Y")
[1] "2012-11-27 PST"

```

So, you need to first convert your date in the appropriate date format and then you will be able to calculate the time interval.

Personally, I prefer to manipulate dates with `lubridate` package because the name of their function is a little bit more intuitive to me. For example, if you have a date to the format "Day/Month/Year", you can use `dmy` function.

With your example, you can do:

library(lubridate)
library(dplyr)

df <- data.frame(Date.Death = c("02/08/2015","02/05/2015"),
                 Date.Operation = c("27/11/2012","27/11/2012"))

df %>% mutate_all(~dmy(.)) %>%
  mutate(Diff_Weeks = interval(Date.Operation,Date.Death) / weeks(1),
                                      Diff_days = Date.Death - Date.Operation)

  Date.Death Date.Operation Diff_Weeks Diff_days
1 2015-08-02     2012-11-27   139.7143  978 days
2 2015-05-02     2012-11-27   126.5714  886 days
```

> 2 votes

---
Tags: r
---