# Question
Title: How should multi allelic sites be handled in calculation per site Fst?
I'd like to calculate per site Fst(fixation index) to investigate genomic differentiation between two populations(haploid).

I did SNP calling and got vcf file. In the vcf file, multi allelic sites exist:

```
#CHROM  POS ID  REF ALT
1   154 .   G   "A, T"

```

For example this site is triallelic sites.

1. Could I calculate Fst in this site by normalizing vcf file?

Normalization means dividing into 2 lows. REF G ALT A(T would be convert G) REF G ALT T(A would be convert G)

In calculation using Popgenome package, we can specify missing sites(Missing sites would not be included in the calculation). In this case, which is better: replacing with 'REF' or replacing with missing sites?

2. Why most studies calculate only biallelic sites?
3. Before analyzing, removing variants based on MAF is make sence? For example, one sites include A(56%), T(44.5%), C(0.5%). In this case, could I remove variant C(replace missing site(NA)) before calculation of Fst?

Best regards,

# Answer
1. Hmmmm .... I personally would not use a single VCF to calculate Fst in the way you are describing. It will also miss one of the nucleotides. I would aggregate across VCFs for a homologous position and use those frequencies. Its designed for heterozygotes, you are describing a haploid genome, like mtDNA (I think). Fst needs a genuine allelic frequency, whereas you've shoehorned it (possibly).
2. Bi-allelic because that is a synonymous variation for most amino acids, i.e. without triggering an amino acid change.
3. In terms of removing a variant ... unless you believe this represents the *in vitro* error rates (which you can define from an invariant site) ... then okay. If it is above the *in vitro* error rate then you can't remove it. I know the read-depth here, so I can't comment.

---

Just to clarify Fst can be used for prokaryotic genomes (haploids), but unless I've misunderstood the question it seems to being derived from a single VCF ... single Metagenome (?). Thats not how it works. I can describe the what the calculation is measuring will artificially inflate Fst and you can't correct for it.
> 0 votes

---
Tags: vcf, snp, bcftools, fst
---


# Question
Title: multiplatt from indicspecies package not producing results for all groups
Based on presence absence data, I'm attempting to find indicator species for different sample sites. I have 25 different groups or "clusters" as I've labelled them. Running the R code below, only results for 13 of the 25 groups was returned.

Is the output truncated, i.e. is there a max output parameter in R I need to fix? Or does this imply that there simply aren't any statistically relevant indicator species in the groups not included in the output?

Code used:

```
library(indicspecies)

data = read.csv("ISAI.txt", header= TRUE, sep= "\t") #Load up data from file

abund = data[,3:ncol(data)] #Collect only numerical data, starts in 3rd column
cluster = data$Cluster #Groups are defined in the Cluster column of the imported data

inv = multipatt(abund, cluster, func = "r.g", control = how(nperm=1000),duleg = TRUE) #Run analysis.

summary(inv) #print out results

```

**Notes**:

For the `func` parameter, "r" was used due to data being in presence/absence format. The "g" corrects for groups with different numbers of members. The `duleg = True` prevents group site combinations from being assessed as I'm not interested in those.

# Answer
> 1 votes
Got it! :)

I switched the `summary(inv)` to `summary(inv, alpha=1)` to get all possible results and all my groups showed up.

So, the answer to my own question is that the groups are missing from the output because they don't have any statistically significant OTUs correlated with those groups.

---
Tags: r, features
---


# Question
Title: How can I filter variants in a VCF file based on size?
I have a filtered VCF file with indels (no SNPs). Does anyone know of a easy way to use bcftools/vcftools to filter out small or large indels (I'm interested in indels between 50-100bp in size)? Thanks!

# Answer
> 1 votes
Vcftools allows you to filter to include or not include a list of sites.

https://vcftools.github.io/man\_latest.html#SITE%20FILTERING%20OPTIONS

If you can first grab the locations of indels of this size (for example you could write a python script that goes through your vcf and grabs the location of anything with ref or alt of a length of at least 50 characters), then you can provide vcftools with a list of variants to filter on.

If you are not comfortable writing your own script to look through the file, you could delete the vcf header lines, change the file extension to tsv, then open it in excel. Then you could use excel to filter for those sites and get their positions.

---
Tags: vcf, genome, bcftools, vcftools
---


# Question
Title: How to extract mutations in a particular gene from the VCF file of whole genome
Also posted on biostars

I would like to extract all the mutations in a single gene of interest from the VCF file. For example I would like to extract all the mutations in the NFKB1 gene from the whole genome VCF file.

# Answer
> 1 votes
Find the location of the gene of interest and then filter the vcf for that location using vcftools, bcftools, plink, or some other software. There is a lot of good documentation available online for these tools.

# Answer
> 1 votes
As mentioned, there are many options to go (eg. vcftools and bcftools). The Broad Institute provides a valuable toolset (GATK) for this very purpose and extended procedures for various variant selection tasks.

SelectVariants function can be good for this purpose, though the location of a gene of interest should be addressed: https://gatk.broadinstitute.org/hc/en-us/articles/13832694334235-SelectVariants

---
Tags: vcf
---


# Question
Title: Feature annotation: RefSeq vs Ensembl vs Gencode, what's the difference?
What are the actual differences between different annotation databases?

My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations. Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

# Answer
> 24 votes
To add to rightskewed answer: While it is true that:

Gencode is an additive set of annotation (the manual one done by Havana and an automated one done by Ensembl),

the annotation (GTF) files are quite similar for a few exceptions involving the X chromosome and Y par and additional remarks in the Gencode file (see more at FAQ - Gencode).

> What are the actual differences between different annotation databases?

They are a few differences, but the main one **for me (and it could be stupid)** is

that **Refseq is developed by the American NCBI** and

the **ENSEMBL is mainly developed by the European EMBL-EBI.**

Often, labs or people will just start using what is the best known to them (because of a course or workshop) or because they start working with one of the databases with one specific tool and keep with it later.

> My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations.

Your lab might be mostly European based people or they might also have read papers like the one from Frankish et al. Comparison of GENCODE and RefSeq gene annotation and the impact of reference geneset on variant effect prediction. BMC Genomics 2015; 16(Suppl 8):S2 - DOI: 10.1186/1471-2164-16-S8-S2

From the Frankish et al. paper paper:

> The GENCODE Comprehensive transcripts contain more exons, have greater genomic coverage and capture many more variants than RefSeq in both genome and exome datasets, while the GENCODE Basic set shows a higher degree of concordance with RefSeq and has fewer unique features.

As for:

> Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

No. I don't think they are great differences between them as that the global picture should stay the same (although you will see different results if you are interested in a small set of genes). However, **they are not directly interchangeable**. Particularly as there are many versions of Ensembl and Refseq based on different genome annotations (and those won't be interchangeable between themselves either in most cases).

**However, you can easily translate** most\[1\] of your Refseq IDs to ENSEMBL IDs and vice-versa with tools as http://www.ensembl.org/biomart/martview for example (there are devoted libraries/API as well like Biocondutor: biomaRt

\[1\] Most as sometimes, they might be annotated in one of the database but haven't (yet) an equivalent in the other.

**EDIT**

In fine, even if people tends to keep to what they are used to (and that the annotations are constantly expanded and corrected) depending on the research subject one might be interested in using one database over another:

From Zhao S, Zhang B. A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics. 2015;16: 97. paper:

> When choosing an annotation database, researchers should keep in mind that no database is perfect and some gene annotations might be inaccurate or entirely wrong. \[..\] Wu et al. \[27\] suggested that when conducting research that emphasizes reproducible and robust gene expression estimates, a less complex genome annotation, such as RefGene, might be preferred. When conducting more exploratory research, a more complex genome annotation, such as Ensembl, should be chosen.
>
> \[..\]
>
> \[27\] Wu P-Y, Phan JH, Wang MD. Assessing the impact of human genome annotation choice on RNA-seq expression estimates. BMC Bioinformatics. 2013;14(Suppl 11):S8. doi: 10.1186/1471-2105-14-S11-S8.

# Answer
> 12 votes
## Ensembl vs Gencode

https://www.gencodegenes.org/faq.html

> The GENCODE annotation is made by merging the Havana manual gene annotation and the Ensembl automated gene annotation. \[...\] In practical terms, the GENCODE annotation is identical to the Ensembl annotation.

Further, for the GTF file differences:

> The only exception is that the genes which are common to the human chromosome X and Y PAR regions can be found twice in the GENCODE GTF, while they are shown only for chromosome X in the Ensembl file.

## Gencode(Ensembl) vs RefSeq

Gencode is in almost all cases more comprehensive. For example, this is NCBI RefSeq vs Ensembl (v24, release 83) for BRCA gene:

RefSeq and Gencode are not interchangeable in most cases, though RefSeq annotations will often be a subset of the Gencode ones.

# Answer
> 6 votes
To add practical advice to what others have said:

In a practical sense, I think the biggest difference between RefSeq and Ensembl/GENCODE is in the sensitivity/specificity trade off.

Ensembl aims more towards the inclusive end, including a far larger number of transcript variants, many of which are only weakly supported.

RefSeq trades some of this sensitivity for specificity - you can be more confident that a RefSeq transcript exists, but less confident that the ReqSeq annotation includes all of the real transcripts for a gene.

# Answer
> 1 votes
While annotations between RefSeq and Gencode are not so different is the coding regions (genes), Gencode is much more rich in the intergenic regions. This could be very advantagous for epigenetic studies, where regulation is of interest.

# Answer
> 1 votes
The UCSC Genome Browser Genes FAQ discusses this question in detail: https://genome.ucsc.edu/FAQ/FAQgenes.html#ens

> Officially, the Ensembl and GENCODE gene models are the same. On the latest human and mouse genome assemblies (hg38 and mm10), the identifiers, transcript sequences, and exon coordinates are almost identical between equivalent Ensembl and GENCODE versions (excluding alternative sequences or fix sequences).
>
> GENCODE uses the UCSC convention of prefixing chromosome names with "chr", e.g. "chr1" and "chrM", but Ensembl calls these "1" or "MT". At the time of writing (Ensembl 89), a few transcripts differ due to conversion issues. In addition, around 160 PAR genes are duplicated in GENCODE but only once in Ensembl. The differences affect fewer than 1% of the transcripts. Apart from gene annotation itself, the links to external databases differ.
>
> The GENCODE Release History shows the release dates and can be linked to corresponding Ensembl releases. You can download the gene transcript models from the website https://gencodegenes.org or from http://ensembl.org. For most applications, the files distributed on the GENCODE website should be easier to use, as the third party database links are easier to parse and the sequence identifiers match the UCSC genome files, at least for the primary chromosomes.
>
> Additional information on this question can be found on the GENCODE FAQ page.

# Answer
> -1 votes
Note that the pDot of a certain mutation might differ based on which transcript you are using, and so Ensembl and RefSeq cannot be used interchangeably.

---
Tags: sequence-annotation, ensembl, refseq, gencode
---


# Question
Title: Gathering data on bacterial organism growth conditions
let's say I have a list of organism names like the example below

```
Achromobacter xylosoxidans
Acidithiobacillus ferrooxidans
Chloroflexus aurantiacus
Clostridium perfringens
Aquaspirillum arcticum

```

The goal is to identify the optimal growth conditions for each organism--meaning the optimal growth temperature and pressure. Obviously I can just google this information if the list is composed of a few organisms or even few tens of organisms. However, assuming that the list is too large, let's say several 100s of names, what would be an efficient means to gathering this information.

In searching for this info, I have not come across a database that has the information that I'm looking for in one place. Instead I find that this information is scattered across multiple sites, databases, scientific publications, protocols, etc. For each organism, I basically have to search through and extract this manually, and I was wondering if anyone can suggest some way to automate this. I am familiar with Python3, ideally I would accomplish this with python.

I though about implementing some sort of a web scraping solution, where I provide the name of the organism, as well as a few key words and see what happens. However, many times this data is not well tabulated and the name and whenever you do find information on the organism, you need to read to be able to interpret the information to get your answer...

maybe I'm not looking in the right place, or my approach to this problem is wrong/backwards/etc, any suggestions will help.

I haven't been able to find a centralized place to gather this information, that's why I posted here. The closest thing I've found is the ATCC.

# Answer
> 1 votes
You could use the BacDive API, or you could download the kgx formatted files from KG Microbe.

---
Tags: database, public-databases
---