cannin/biostars_qa huggingface dataset:
Question: I might misunderstanding the following idea:

Starting at the variant of interest, the region is extended out ±0.1cM, using HapMap fine-scale estimates of recombination rates.

From the HapMap you get file containing the following:

```
chrom start stop Rate_cM.Mb Avg_cM Gen_map_cM
Chr1 45413 72433 2.27700294070856 0.061526896460886 0
Chr1 72434 78031 2.42987467531016 0.0136024384323863 0.061526896460886
Chr1 78032 227743 2.43386172133186 0.364378306024035 0.0751293348932723
```

If a variant falls into 45413 - 72433 range then in order to calculate +- 0.1 cM region around it, I need to proceed as follows:

72411-45413 = 27009 bp , this region has an average 0.0615 cM, meaning that 0.1cM corresponds to 43917 bp ( 43917 = 27009*0.1/0.0615). Is it right?

Thank you in advance.

Answer: Not quite. Remember that the recombination rate changes with position, so the expected recombination rate over the next 43917 bases isn't going to be 0.1. So if we're interested in +0.1cM from 72411, then getting to 72433 gets us 4.844534e-05cM. Getting to 78031 adds an additional 1.36e-2, bringing us to 0.01365088. The means we need another 0.08634912cM. The last entry has an average rage of 2.434e-6cM/base, so dividing gives us an additional ~35478 bases. That puts us at 113509, meaning that the total width is 41098 bases.

The distance also won't be symmetric (i.e., -0.1cM probably won't be 41098 bases in the other direction).
Question: Hello folks,

I would like to use a gene name (*e.g. DOCK2*) or its ENSEMBL ID as input, and get as output the reported GWAS traits to the respective gene. I did give a look into the [GWAS Catalog API](https://www.ebi.ac.uk/gwas/docs/api) webpage, but I didn't find a way to use gene as input. Also, I looked for R and Python libraries for doing that, and I couldn't find either (maybe I didn't search well).

I know it's possible to do that for associations, studies, SNPs, and efoTraits. For instance, I could just use https://www.ebi.ac.uk/gwas/rest/api/singleNucleotidePolymorphisms/rs4918943 url, where it returns info for the SNP "rs4918943". I would like something similar for a gene ID, or a library that does it.

It would be great to input a gene and receives any output where I could retrieve the GWAS trais related to such gene.
Thanks in advance!

Answer: I didn't find the usage for it, but I figured out how to query it by using the following URL:

```
https://www.ebi.ac.uk/gwas/api/search/downloads?q=ensemblMappedGenes:{gene}&pvalfilter=&orfilter=&betafilter=&datefilter=&genomicfilter=&genotypingfilter[]=&traitfilter[]=&dateaddedfilter=&facet=association&efo=true
```

where we should replace `{gene}` with the gene symbol (*e.g.* TP53).

**Example using R:**
```r
# Required library # could be readr::read_delim() as well
library(data.table)

# Set the function
genesymbol2gwas <- function(gene){
url <- paste0(
"https://www.ebi.ac.uk/gwas/api/search/downloads?q=ensemblMappedGenes:", gene,
"&pvalfilter=&orfilter=&betafilter=&datefilter=&genomicfilter=&genotypingfilter[]=&traitfilter[]=&dateaddedfilter=&facet=association&efo=true"
)
return(fread(url))
}

# Use the function for gene "TP53"
genesymbol2gwas("TP53")
```

**Example using Python:**
```python
# Load library to read the result table as a dataframe
import pandas as pd

# Define the function
def genesymbol2gwas(gene):
url = "https://www.ebi.ac.uk/gwas/api/search/downloads?q=ensemblMappedGenes:{}&pvalfilter=&orfilter=&betafilter=&datefilter=&genomicfilter=&genotypingfilter[]=&traitfilter[]=&dateaddedfilter=&facet=association&efo=true"
return pd.read_csv(url.format(gene), sep='\t')

# Use the function for gene "TP53"
genesymbol2gwas("TP53")
```

Hope it helps someone in a similar situation =)
Question: Hi all,

I am using shapeit to prephase the genotypes, followed by using minimac to do the imputations (with 1000g phase1 v3 as the reference panel)

Now when I phase the genotypes - I would like to use the genetic map from 1000g instead of the hapmap genetic map that is mentioned in shapeit manual..

Is it possible for someone to direct me to the right source.. also it would be great if you can correct if my reasoning is wrong.

Thanks

Answer: Have you checked the Impute2 1000 Genomes phase 3 reference data?

https://mathgen.stats.ox.ac.uk/impute/1000GP_Phase3.html
Question: I wanna to see the expression profile in the ucsc genome browser of several genes. I was thinking if someone know some way or script that can retrieve the png images automatically from the genome browser.

Answer: I posted a Python-based tool on Github called [soda.py][1] that creates a web-ready gallery of UCSC browser shots. You just give it your BED file of coordinates, build, and session ID, and you specify an output directory where PDF and PNG results get stored (as well as an `index.html` file that lets you browse through snapshots with a web browser).

If you want to do things by hand, you can do something like the following quick and dirty approach to get a nice PNG. You'll need ImageMagick `convert` installed in order to convert the PDF to PNG. You'll also need GNU `wget` to do web requests on the command line.

$!/bin/bash

chrom="chr1"
chromStart=1234567
chromEnd=1234987
sessionID=1234
genomeBrowserURL="genome.ucsc.edu"
dumpURL="http://${genomeBrowserURL}/cgi-bin/cartDump"
postData="hgsid=${sessionID}&hgt.psOutput=on&cartDump.varName=position&cartDump.newValue=${chrom}%3A${chromStart}-${chromEnd}&submit=submit"
wgetOpts="--no-directories --recursive --convert-links -l 1 -A hgt_*.pdf"
wgetWaitOpts="--wait=1 --random-wait --tries=2 --timeout=100"
wget ${wgetWaitOpts} --post-data "${postData}" "${dumpURL}"
wget ${wgetOpts} ${wgetWaitOpts} "$url&position=${chrom}%3A${chromStart}-${chromEnd}" 2> fetch.log
mv hgt_*.pdf ${sessionID}.pdf
convert -density 300 ${sessionID}.pdf -background white -flatten ${sessionID}.png

You'd fill out `chrom`, `chromStart`, `chromEnd` and `sessionID`. Or use placeholders `$1` etc. and pass them in on the command line.

If you have a few regions to look at, this would just be a matter of modifying this approach to loop over their respective chromosome name and interval values, and naming the output files appropriately. (Or you can use [soda.py][1] for automation.)

[1]: https://github.com/alexpreynolds/soda
Question: Hi, Im a begginer in MEGA software. I have 89 protein sequence for which I need to construct a phylogenetic tree using bootstrap method with 1000 replication with data set parameter with complete deletion. But I am not able to construct a tree because of 3 sequence whose protein length is very less when compared to other 86 sequence. Even I tried by deleting non conserved regions in all protein sequence but still I am not able to get a tree because the size of the smaller proteins become smaller and smaller. Kindly help me out in solving this problem.

Answer: Regardless of the approach or program you are using, the input for any phylogenetic estimation approach is an alignment, i.e., an inference of homology. Therefore, by necessity, your sequences **must have a shared ancestry** to even begin to infer a phylogeny. If the sequences are shorter but homologous, a multiple sequence alignment (of nucleotides or amino acids or both via a translation alignment for protein-coding sequences) ought to resolve the sequences by introducing gaps - insertions or deletions. It sounds like you're not doing this; when you say

> The 3 short protein sequence are upregulated in abiotic stresses. Is it ok if I omit the sequence because they have role in abiotic stresses?

It suggests that your dataset may consist of multiple proteins, not the same protein across samples, which is a **completely inappropriate** input for phylogenetic techniques.

In other words, your workflow would be:

1. Construct a dataset of the same locus across all samples
2. Align the amino acids or nucleotides
3. Model selection for ML analysis or NJ distance corrections/uncorrected NJ/UPGMA/etc.
4. [If you decide to use a model: With an appropriate model, any likelihood (maximum likelihood or Bayesian) approach.]
5. Bootstrapping etc. for support.

If you do have sequences with a shared history, I would follow Istvan Albert's recommendation and remove the short sequences if they are truly unalignable.
Question: Probably a basic question but it's difficult to find information on the subject...

When viewing a GTF and BED file in IGV, there seems to be differences in colour (BED track can't be recoloured, just shows as black, GTF can be recoloured). Does anyone know the reason for this? I thought that the GTF and BED formats contained predominantly the same information just in different formats.

Thanks!

Answer: BED was originally designed to contain color specification for each row in column 9. In addition it can contain a number of track definitions in the header that can further specify the colors. Using these have fallen out of favor (thankfully so I might add) but IGV seems to stick with the more strict definition of reading the color from the rows themselves.
Question: In the last few weeks I have been using [agriGO][1] to perform GO enrichment analysis on a non-model organism dataset. In the last week or so I've noticed that the website hosting it is down. I know it is not just my connection since https://isitup.org/bioinfo.cau.edu.cn also shows that is is down. Does anyone know what the status of agriGO is? If it will come back online or if there is a mirror?

Alternatively, I'm looking for a similar service that I could supply a background list of GO terms, and a regulated list of GO terms to perform enrichment analysis of.

[1]: http://bioinfo.cau.edu.cn/agriGO

Answer: I guess they turned off agriGO v1 (yours, paper from 2010) and left only agriGO v2 (paper from 2017) up? http://systemsbiology.cau.edu.cn/agriGOv2/

Especially for GO annotation I would always use the newest data possible - see for example https://www.nature.com/articles/s41598-018-23395-2
Question: Hello, I am trying to find where to download human alpha satellite sequence. So far I have only been able to find one from 1987 paper: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC306152/pdf/nar00260-0438.pdf

I know that Karen H. Miga incorporated her Centromere Reference Models into GRCh38, so downloading her predicted sequences of centromeres would also be helpful - if you can point me where to download those. However, alpha satellite (or multiple since they are variable) would be the best.

Answer: Have you tried [RepBase](http://www.girinst.org/server/RepBase/index.php)?
Question: Hi friends,

I am working on a population genetic project using microsatellite markers. Due to the large sample size and 12 microsatellite loci, I will need to use an efficient and solid method/tool to score the fragment analysis data (.fsa file type) from genetic analyzer 3130. In total, I will have nearly 5000 files to be scored. I have done part of my pilot fragment analysis results with a free version software, PeakScanner 1.0 (ABI); however, this software doesn't come with automatic features, and thus I have to score the data manually, which is very time consuming and inconvenient. In addition, we don't have the budget to purchase some commercial softwares (such as, Gene Mapper ) to do this. I am wondering if you guys can give me some suggestions to handle this situation.

OR

I can keep the question short, any substitution software package for ABI GeneMapper to score multiplex fragment analysis data (*.fsa)?

Thank you very much.

Answer: Update for this question:

There's a package called `Fragman` in R just released in 2015, performs fragment analysis using FSA files: https://cran.r-project.org/web/packages/Fragman/
Question: Hi,

I have SNP data for approximately 600,000 SNPs that I'll be using for eQTL analysis.

I've been advised to derive a pruned set of SNPs that are in approximate linkage disequilibrium (LD).

I've used the SNPRelate::snpgdsLDpruning function in R to do this, using an LD threshold of 0.2.

Is this an appropriate threshold to use?

This takes our 600,000 snps down to ~60,000; using a threshold of 0.1 leaves ~20,000 snps.

This is quite a big difference so I'm wondering if 0.1-0.2 is too stringent.

Is there a more-or-less standard threshold that is used for LD pruning?

Thank you.

Answer: The choice of optimal r^2 threshold for LD pruning is highly dependent on the population history of your study subjects. Despite D' that is purely a measure of the non-random association of alleles at two or more loci, r^2 values are informed by allele frequencies as well. For example for two polymorphic loci (A & B), one with 50% allele frequency and the other with 1% frequency which are in complete LD, the D' value would be 1, but r^2 would only be 0.01. This tells us although these two polymorphic loci are in complete LD with each other, the allele B is so rare that 99% of the time it is not observed on the same haplotype with allele A.

To circumvent such situations, a minor allele frequency filtering is carried out prior pairwise LD calculation to remove rare alleles. Sometimes, if your study population is not a true representative of the extant population, this filtering step is not sufficient since rare alleles can rapidly drift up to higher frequencies in structured or founder populations. Although I assume these scenarios are an overstretch to your question, you'd like to use more stringent r^2 values to avoid colinearity of effects among your pruned SNPs.

On the contrary, if you're dealing with a population with massive haplotype diversity (such as sub-Saharan Africans) and you would assume that your study population is not large enough to be a true representative of all haplotypes in the population you may want to use a more relaxed r^2 threshold for pruning.

Overall, as long as you can rationalise your choice of r^2 either thresholds are acceptable, but r^2 < 0.2 is the common practice for European populations.
Question: I am trying to annotate a variant file(generated using strelka) from mice WGS data. This is the command I used:

./vep -i /path/to/somatic.snvs.vcf \
--cache /data/shayantan/mus_musculus/ \
--species mus_musculus

The output variant file has no gene names. Why is this happening? Something wrong with my cache files?

EDIT (@Ram): Sample input VCF:

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NORMAL	TUMOR
chr1	3003110	.	G	T	.	LowEVS	SOMATIC;QSS=17;TQSS=1;NT=ref;QSS_NT=17;TQSS_NT=1;SGT=GG->GT;DP=35;MQ=60.00;MQ0=0;ReadPosRankSum=1.95;SNVSB=3.58;SomaticEVS=0.80	DP:FDP:SDP:SUBDP:AU:CU:GU:TU	18:0:0:0:0,0:0,0:17,17:1,1	14:0:0:0:0,0:0,0:11,13:3,4
chr1	3035137	.	G	A	.	LowEVS	SOMATIC;QSS=17;TQSS=2;NT=ref;QSS_NT=17;TQSS_NT=2;SGT=GG->AG;DP=70;MQ=40.40;MQ0=10;ReadPosRankSum=0.89;SNVSB=3.23;SomaticEVS=0.10	DP:FDP:SDP:SUBDP:AU:CU:GU:TU	27:0:0:0:3,7:0,0:24,28:0,0	27:0:0:0:4,6:0,0:23,29:0,0
chr1	3035168	.	C	T	.	LowEVS	SOMATIC;QSS=8;TQSS=2;NT=ref;QSS_NT=8;TQSS_NT=2;SGT=CC->CT;DP=51;MQ=47.72;MQ0=3;ReadPosRankSum=1.78;SNVSB=2.68;SomaticEVS=0.08	DP:FDP:SDP:SUBDP:AU:CU:GU:TU	18:0:0:0:0,0:16,19:0,0:2,4	23:0:0:0:0,0:20,25:0,0:3,3
chr1	3035504	.	C	A	.	LowEVS	SOMATIC;QSS=15;TQSS=2;NT=ref;QSS_NT=14;TQSS_NT=2;SGT=CC->AC;DP=59;MQ=51.03;MQ0=2;ReadPosRankSum=-1.19;SNVSB=2.71;SomaticEVS=0.09	DP:FDP:SDP:SUBDP:AU:CU:GU:TU	23:0:0:0:3,5:20,22:0,0:0,0	27:0:0:0:4,4:23,28:0,0:0,0
chr1	3043000	.	G	T	.	LowEVS	SOMATIC;QSS=21;TQSS=1;NT=ref;QSS_NT=21;TQSS_NT=1;SGT=GG->GT;DP=53;MQ=46.60;MQ0=7;ReadPosRankSum=1.70;SNVSB=1.37;SomaticEVS=0.20	DP:FDP:SDP:SUBDP:AU:CU:GU:TU	20:0:0:0:0,0:0,0:18,24:2,3	22:0:0:0:0,0:0,0:18,22:4,4

Answer: [Those variants are all intergenic][1]. There is no gene symbol because no genes are hit.

EDIT (@genomax) - Actual answer is further below in this chain at https://www.biostars.org/p/331209/#331504

[1]: http://www.ensembl.org/Mus_musculus/Tools/VEP/Ticket?tl=F9VQxvXxHkzscTq9
Question: Hi

Can we sort BAM files according to the read name? As the normal sorting happens on the co-ordinates, can anyone tell how to sort a BAM file on read names? I am trying to run HT-Seq count on paired end SAM files but receiving warnings for which I have to sort the BAM in read names and then create its SAM and then run HT-Seq

Answer: See Picard's SortSam command below. You can sort a sam or bam file using queryname or readname.

http://picard.sourceforge.net/command-line-overview.shtml#SortSam
Question: Hello,

I used the RidgePlot to make the plot. What is the definition of "Identity", what is the mean of the "identity"? Is that possible I can delete "Identity"?

RidgePlot(Plasma, features = c("UCP2"), group.by="tech", ncol = 3)

Thank you in advance for your great help!

Please see attached file

Thank you!

Best,

Yue

Answer: Hi,

`Identity` is the identity of the cells in scRNA-seq data. It's not directly related with this plot, but with scRNA-seq data in general. So, the `identity` is any factor/categorical variable that describes each cell. The `identity` of the cells can be the sample, replicate, condition, cell cluster to which they belong to.

In your case it seems to be `condition` or `sample`. But you can set the `identity` to cell clusters/populations/types any variable that annotates each cell to a factor/categorical variable.

You can see the `identity` of your cells by typing:

Idents(pbmc) # assuming that pbmc is your seurat object

You can set the `identity` of your cells to another variable in your seurat object by doing:

Idents(pbmc) <- "replicate"

Where `replicate` is a variable in your seurat object - "pbmc" - annotating each cell as `rep1` or `rep2` (see the example here: https://satijalab.org/seurat/v3.0/interaction_vignette.html ).

I hope this answers your question,

António
Question: Hi all,

I have RNA-seq data for 3 conditions (control, breast, endometrial). Some of the samples where collected and sequenced in two different batches. I want to correct for batch effects using EdgeR and GLMs. I did read the EdgeR vignette, and specifically section 4.5 but I'm still a bit confused with the contrasts matrix and if I'm interpreting the results correctly.

This is my code:

design = model.matrix(~group+batch, data=d$samples)
d$samples$group = relevel(d$samples$group, ref="control")
rownames(design) = colnames(d)

d <- estimateGLMCommonDisp(d, design, verbose=TRUE)
d <- estimateGLMTrendedDisp(d, design)
d <- estimateGLMTagwiseDisp(d, design)

fit = glmFit(d, design)
lrt = glmLRT(fit, contrast=c(0,0,-1,0))

This is my design matrix

intercept	groupbrc	groupendo	batch2
sample.1	1	0 1 0
sample.2 1 0 1 0
sample.3	1	0 1 0
sample.4	1	0 0 0
sample.5	1	1 0 0
sample.6	1	1 0 0
sample.7	1	1 0 1
sample.8	1	1 0 1

Basically what i want to do is pairwise comparisons between the treatments (brc vs normal, endo vs normal, brc vs endo) and accounting for batch effects at the same time. I understand that the (intercept) corresponds to the normal condition but what I don't understand is what the last column (batch 2) means, and if i should include it in my contrasts. The contrasts I've used are the following, i get DE genes but I'm not sure if I'm accounting for batch effects correctly

Brc VS normal `lrt = glmLRT(fit, contrast=c(0,1,0,0))`

Endo VS Normal `lrt = glmLRT(fit, contrast=c(0,0,1,0))`

Brc VS Endo `lrt = glmLRT(fit, contrast=c(0,1,-1,0))`

Thank you!

Answer: Brc vs. normal: `lrt = glmLRT(fit, coef=2)`
Endo vs. normal: `lrt = glmLRT(fit, coef=3)`

I would do **Brc vs. Endo** the same as you did, since you do want to contrast two model coefficients. There's no need to use a contrast when you're just wanting to test whether a model coefficient is itself significant.
Question: Hello!

I am very new to biopython and I am trying to accomplish what I think is a simple task: I would like to remove sequences from a protein alignment that do not contain a particular residue at a specified position. I would like to be able to input a protein alignment in fasta format and then output a new alignment where all the sequences that do not meet my criteria are removed

For example: My input protein alignment contains sequences that have a mixture of residues at position 137. I would like to output a new alignment that contains only sequences that have either an arginine or a valine at position 137.

Just a bit of additional clarification: I am sequencing an amplicon of a functional gene and generating protein sequence alignments using RDP's fungene pipeline. I want to further screen the alignment by eliminating any sequences that do not contain a selection of conserved residues at various positions.

Thank you very much for your time.

-J

Answer: python script.py

- read in the sequence `from Bio import SeqIO` and store them in a dictionary. Print Sequences in a file.
- align the sequence file with tool of your choice.
- read in the alignment using from Bio import AlignIO
- iterate through the alignment, and check the residues you are interested.
- make a list of the ones to kept or thrown. Delete your alignment file.
- Print the remaining sequences and re-do alignment.

Hope I haven't confused it.
Question: Hi, all!

I'm trying to do the pathway analysis from bovine RNA-seq data using 'gage' in R, but I couldn't find the 'kegg gene set' for bovine.

If anyone knows how to obtain it, please help me!

Thank you!

Answer: In gage package, you can use kegg.gsets function to generate kegg pathway gene set data or ~3000 KEGG species, and go.gsets function to generate GO gene set data for 19 major species.
Unlike the ones provides with gage and gageData, these gene set data are most up-to-date ones.
For details, check:

?kegg.gsets

?go.gsets

The main vignette of gage also showed some example in these functions:
http://bioconductor.org/packages/release/bioc/vignettes/gage/inst/doc/gage.pdf
Question: What are ‘reads’ in mean reads per cell? What are 'reads' in mean reads per cell?

Answer: If this question is about 10x `cellranger` data then details can be found [at this link][1].

Mean Reads per Cell	= The total number of sequenced reads divided by the estimated number of cells.

[1]: https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/metrics
Question: Hi all,

I know this may sound simple, but I am really having trouble finding a way how to extract reads that fall **outside** a give region. I tried the -U option but it doesnt work (or maybe I am not using it correctly). I didn't find any example either which uses -U option of samtools view.

I am using this command:
$samtools view sample.bam 2:33050509-33154206 -U without-region.sam

This is paired end data and I want to retain the other paired read that falls outside the region, therefore grepping reads other than the one in my specified region, wont work either.

Thanks in advance.

Answer: Hi all,

I found the answer to this issue, its a very simple fix. Sorry for replying so late. Got busy with some other stuff and totally forgot to reply. Anyways, if one wants to select regions outside of a given region, you can use -U option, but as I said earlier it doesnt work with the range specified on command line. But if you give a bedfile containing the region ad then provide -U option, it will work.

For example:
Create bedfile (my.bed):
2 33050509 33154206 LINC00486
and now use this bedfile to extract region other than the ones in the given bedfile by:
$samtools view sample.bam -L my.bed -U without-bed > /dev/null

To my knowledge with this approach, a read that has even a single base in the specified region gets counted as within region, the read just need to span over either of the start or stop of the region (even if its a single base).

Let me know if you have any questions.

Thanks.