cannin/biostars_qa huggingface dataset:
Question: When I'm working with Mouse ChIP-seq data, I normally remove mapped reads which overlap the ENCODE blacklist regions. Previously there was no data for the mm10 assembly, so instead I would lift over the coordinates from the mm9 assembly ([as described in the F1000 csaw article][1]). When I do this I get 3,010 regions. However, recently ENCODE created a dataset for the mm10 assembly, but it only contains 164 regions. I contacted ENCODE to ask why and this was their response:

"LiftOver is not a good strategy for transferring blacklists across assemblies. Note that the blacklists are regions that show artifacts due to deficiencies in the genome assembly (e.g. unannotated repeats). So with a better assembly a region that was previously a blacklist wont be one any more. GRCh38 and mm10 have fewer detectable artifacts compared to GRCh37 and mm9 respectively because they are better and more complete assemblies e.g. repeats near centromeres and telomeres are better annotated. Hence the fewer regions. This blacklist release is also a first pass for mm10 and GRCh38 with minimal manual curation. We will be releasing additional refined versions in the future that may capture additional regions."

This made me doubt the advice given in the csaw paper, and my usual processing stages. I've always followed the advice from Heng Li that you should [map to an un-masked genome][2]. Then I usually remove the liftOver ENCODE blacklist regions. Instead should I use ENCODE's official mm10 blacklist, and then also remove predicted repeat regions from the UCSC genome annotation?

[1]: https://f1000research.com/articles/4-1080/v2
[2]: https://www.biostars.org/p/3232/

Answer: We've usually taken the approach of just removing peaks near blacklisted regions (we've been using a list from the DKFZ with ~400 regions) when doing peak calling and otherwise just instructing tools to ignore those regions when that's possible (e.g., in deepTools).

As an aside, I agree with the ENCODE response regarding the problematic nature of lifting over blacklisted regions. As an example, only 10% of the problematic regions in hg18 were still problematic in hg19. Somewhere Heng Li has a presentation showing further improvements in GRCh38. I would presume this is the same for the mm8 -> mm9 -> mm10 progression of releases.
Question: Hi friends,

I am trying to download ht-seq file form GDC portal but it gives me one file for each patients. Can you please let me know how to download one file including all patients together for all 60000 genes?

Is there any R code for this?

Thank you

Answer: You can use:

[xenabrowser][1]

or

r package [TCGAbiolinks][2]

[1]: https://xenabrowser.net/datapages/
[2]: https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html
Question: We are trying to integrate STAR into a galaxy pipeline. As you might know, STAR can only run as a single process per directory, which is currently the main obstacle. The main reason for star not being able to run in parallel is that names of the output files are not parameters but fixed. The program will refuse to run even when detecting an already running star in the directory.

Is there a plan or a way to make the output file names configurable using a prefix so it will not fail? Maybe a source modified version with this feature?

Second best option, to give each star its own chrooted tmp environment, but I would like to go for the cleaner option first.

Answer: I've gotten around this on our my old cluster by giving each instance a different --outFileNamePrefix, does that not suffice for you? This does still create some annoying empty temp files in whatever directory you execute STAR in that need to then be removed (they don't block multiple instances though).
Question: Hi I am trying to annotate the variants in a vcf file using snpEFF, but its giving multiple annotations for a single position. Any suggestions on how to choose the best or most appropriate one?

Thanks in advance

Answer: Firstly, I'd suggest checking some of the results with another annotation tool. Perhaps get a few of those variants giving multiple annotations and confirm the pattern with the rather popular and powerful [Variant Effect Predictor][1], or VEP. As suggested already in the comments here, the multiple annotation will be due to alternatively splicing and there is nothing much we can do apart from trying to chose one or a few "best" ones based on certain parameters and based on the question you want to answer, really. Let's suppose you are interested in the most severe consequence of a given variant. If so, choose the most severe consequence with the VEP(--[most_severe][7] flag). If you want to focus on the reliability of a transcript, you can pick the annotations for [Ensembl genes][2] that are golden (identical annotation between Ensembl and HAVANA). You can also filter out those that have lower [TSL][3] (use the [--tsl][4] flag on the VEP script) or choose the 'canonical' (find out what is deemed as [canonical transcript][5] by Ensembl) transcript (use the [--canonical][6] flag), and so on so forth...

[1]: http://www.ensembl.org/info/docs/tools/vep/index.html
[2]: http://www.ensembl.org/Help/Glossary?id=78
[3]: http://www.ensembl.org/Help/Glossary?id=492
[4]: http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html#opt_tsl
[5]: http://uswest.ensembl.org/Help/Glossary?id=346
[6]: http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html#opt_canonical
[7]: http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html#opt_most_severe
Question: I'm trying to convert Ensembl IDs to Gene symbols within a summarized experiment object (more or less an expression matrix) using BioMart.

mart <- useDataset("hsapiens_gene_ensembl", useMart("ENSEMBL_MART_ENSEMBL"))
genes <- rownames(gse_cellgenefiltered_cohort1)
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id", "hgnc_symbol"),values=genes,mart= mart)

For some reason, there is a discrepancy between the number of Ensembl IDs I supply BioMart with and the number of Ensembl IDs it returns.

length(rownames(gse_cellgenefiltered_cohort1))

>[1] 23395

length(G_list$ensembl_gene_id)

>[1] 23316

Another thing I noticed, is that BioMart returns duplicated Ensembl IDs for some of them.

length(unique(G_list$ensembl_gene_id))

>[1] 23314

I don't think there are any duplicated Ensembl IDs in the expression matrix.

length(unique(rownames(gse_cellgenefiltered_cohort1)))

>[1] 23395

Would anyone know why this might be happening?

Answer: There's a few things that might be going on, and it's hard to tell exactly without some examples of the missing or duplicated gene IDs. Here's some ideas though.

BioMart will silently drop any element in `values` that aren't found in the query. There's no error or anything, you just don't get a hit. That's easy to see with a single value, harder to spot in 23,000:

```
## query not found in Ensembl
getBM(values = c("ENSG_NOT_REAL"),
filter = "ensembl_gene_id",
attributes = c("ensembl_gene_id", "hgnc_symbol"),
mart = mart)
#> [1] ensembl_gene_id hgnc_symbol
#> <0 rows> (or 0-length row.names)
```

You can try to identify what input values aren't returned in the results with something like `genes[ !genes %in% G_list$ensembl_gene_id ]`. If that finds something I'd search the Ensembl website manually with a few of the IDs and try to understand why they might be missing from BioMart e.g. they might be from an old Ensembl version and have been retired - there are probably many possible reasons.

For completeness I'll also point out that Ensembl BioMart will ignore duplicate entries in the the `values` argument e.g..

```
## duplicated input values
getBM(values = c("ENSG00000010404", "ENSG00000010404"),
filter = "ensembl_gene_id",
attributes = c("ensembl_gene_id", "hgnc_symbol"),
mart = mart)
#> ensembl_gene_id hgnc_symbol
#> 1 ENSG00000010404 IDS
```

However it looks like you've already checked this isn't the case in your data.

Regarding the duplicated entries in the results, this can occur if there is a one-to-many mapping between the two ID types you're trying to find e.g.

```
## one-to-many mapping
getBM(values = "ENSG00000277796",
filter = "ensembl_gene_id",
attributes = c("ensembl_gene_id", "hgnc_symbol"),
mart = mart)
#> ensembl_gene_id hgnc_symbol
#> 1 ENSG00000277796 CCL3L3
#> 2 ENSG00000277796 CCL3L1
```

Mapping between IDs from different organisations is never perfect and it's pretty common to see instances like this, where a single Ensembl ID maps to two HGNC IDs (or vice versa). You could try to identify the duplicated entries with

```
G_list[ duplicated(G_list$ensembl_gene_id) | duplicated(G_list$ensembl_gene_id, fromLast = TRUE), ]
```
Question: Hi all, I have recently started a project on modelling a polygenic risk score model to evaluate its utilitiy in predicting a certain disease. After doing some reading, I have come across various models for unweighted and weighted Genetic Risk Score models.

I am wondering about the use of odds ratio (OR) versus the Beta-coefficient of each SNP variant in a risk score model. For instance, [here][1] they used the Beta-coefficient in their model, while [here][2] they used the odds ratio. Is there any difference in using the odds ratio versus the Beta-coefficient in a risk score model? Also, I noticed that some papers use log(OR) rather than ln(OR), is there a major difference between both?

Thanks!

[1]: https://journals.plos.org/plosmedicine/article?id=10.1371/journal.pmed.1002383
[2]: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0131603

Answer: Hey,

The odds ratio (OR) is the exponent of the beta coefficient. The beta coefficient itself is the per unit increase/decrease in the exposure. A practical example will explain it better:

# create random data with a condition A + B ('outcome') and *gene1* ('exposure')

modeling <- data.frame(
condition=factor(c(rep("A",100), rep("B",100)), levels=c("A", "B")),
gene1=c(runif(100), runif(100)))
head(modeling)
condition gene1
1 A 0.3607443
2 A 0.3268301
3 A 0.4237005
4 A 0.7621534
5 A 0.1456797
6 A 0.3201094

Note that we have set `A` as the reference level.

# create a binomial logistic regression model, with *gene1*'s expression 'predicting' the outcome

model <- glm(condition ~ gene1, data=modeling, family=binomial(link='logit'))
summary(model)

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.03767 0.30014 -0.126 0.900
gene1 0.07294 0.51258 0.142 0.887

Here, the beta coefficient for *gene1* in relation to condition B versus A is 0.07294. So, *gene1* increases in expression in condition B (if it decreased, the beta coefficient would be negative). This is not a statistically significant finding, though, with *p*=0.887.

We can also test the *gene1* via the Wald test on the beta coefficient:

require(aod)
wald.test(b=coef(model), Sigma=vcov(model), Terms=2)
Wald test:
----------
Chi-squared test:
X2 = 0.02, df = 1, P(> X2) = 0.89

# obtain the odds ratio and upper / lower confidence intervals (CIs):

exp(cbind(OR=coef(model), confint(model, level = 0.95)))
OR 2.5 % 97.5 %
(Intercept) 0.9630288 0.5333548 1.737047
gene1 1.0756700 0.3930383 2.949061

So, odds ratio is just 1.1, which, as you can tell, is not huge and only reflects a slight increase.

# log OR
The log OR is just the natural logarithm of the OR. With regard to why we may even want to use log OR over OR, well, there are probably many reasons. One is that we can calculate the Z score from the log OR:

OR <- 1.0756700
lowerCI <- 0.3930383
upperCI <- 2.949061

logOR <- log(OR)
logORlowerCI <- log(lowerCI)
logORSE <- (logOR - logORlowerCI) / 1.96

Then calculate Z:

logOR / logORSE
[1] 0.1420052

Kevin
Question: My aim is to perform **taxonomic classification and transcript abundance estimation** of **viral** communities from metagenomics sequenced reads. As the first step, I mapped all the reads to a reference which was built (bowtie2-built) based on refseq [viral genomes][1]. Currently, I have SAM and BAM files, and I am looking for tool(s) which can perform the profiling task.

I have came across [viRome package ][2], but would like to know whether more tools/methods are available.

[1]: https://www.biostars.org/p/13646/
[2]: http://www.ark-genomics.org/bioinformatics/virome

Answer: If you are looking for a k-mer based tool then checkout: [DisCVR][1] and if you are interested in building de novo assemblies then you can use: [MetaViC][2]
You might find this paper useful: [Bioinformatics tools for analysing viral genomic data][3]

[1]: http://bioinformatics.cvr.ac.uk/discvr.php
[2]: https://github.com/sejmodha/MetaViC
[3]: http://boutique.oie.int/extrait/20orton271285.pdf
Question: I couldn't find a proper documentation of the software used for generating the read counts of the TCGA level 3 data.

I have done 21 normal sample Vs 21 tumor samples analysis using TCGA RNASeq level3 data to find deferentially expressed genes using DESeq.

And further I have taken Illumina body map SRA file, processed using TOPHAT and generated counts using HTSeq. The HTSeq read counts generated using TOPHAT bam files were compared with 21 tumor sample from TCGA level 3 data.

So, now as expected the differential expressed genes using DESeq between "Illumina body map comparison with 21 tumor samples" and "21 normal sample Vs 21 tumor samples from TCGA" should have good overlap of deferentially expressed genes. But the overlapping genes are very less.

Does this means there is something wrong with the processing of illumina body map file and or due to the variation in protocol followed for TCGA data?

Could anyone tell me how the read counts in the TCGA level 3 data is generated? Using which program?

Answer: TCGA V1 analysis (old) used BWA and the V2 analysis (new) which uses MapSplice. All V1 data was reprocessed as V2.

[This file][1] has details about the V1 and V2 pipelines

Note: Since TCGA data has now moved to GDC, the description file has been updated to point to the new resource.

From the description file

```
V1_BWAtoTranscriptome, V1_RNASeqQuantification: UNC V1 RNA-Seq Workflow - BWA Alignment to Transcriptome
Date: 20101108
```

And

```
V2_MapSpliceRSEM: UNC V2 RNA-Seq Workflow - MapSplice genome alignment and RSEM estimation of GAF 2.1
Date: 05-10-2012
```

[1]: https://gdc-api.nci.nih.gov/legacy/data/845ae614-d24e-4fa4-b1a8-05e0845f288d
Question: Hi guys, is there a way to plot the average expression of a list of genes on the y axis of a volcano plot instead of -log10(FDR/P-val) using EnhancedVolcano or another way (like ggplot2)?

Thank you in advance

Answer: Yes, this is then called an MA-plot, where usually the average expression is on the x-axis and the logFC is on the y-axis.

See https://www.biostars.org/p/457769/
Question: Hello,

We search a kind of score whilch can eliminates many FP call in CNVkit results. We think that the weight can be used.
We search a signification of this weight, an equation. In the manual of CNVkit, we found this extract :

A weight is assigned to each remaining bin depending on:

1. The size of the bin;
2. The deviation of the bin’s log2 value in the reference from 0;
3. The “spread” of the bin in the reference.

(The latter two only apply if at least one normal/control sample was used to build the reference.)

So, we have many questions :

- what "bin" means ? The segment ? The bin in my_target.bed ? my_antitarget.bed ?
- what is the "spread" ? the number of bins in the segment ? the length of the segment ? something else ?
- the score seems to be strongly affected by segment size and / or bin number in the segment, is this the case?
- is it possible to get an equation of the weight ?

Thank you

Alice

Answer: So:

- Bins are the unsegmented regions seen in my_target.bed, my_antitarget.bed, and the .cnn and .cnr [files emitted by CNVkit][1]. Sometimes also called "probes" in the code.
- Spread is the statistical spread of coverages in a bin observed across all of the samples in your [pooled reference][2], similar to standard deviation but calculated differently to be more robust to outliers.
- The `weight` listed in the segmented .cns files is the sum of the weights of the bins/probes spanned by the segment. It correlates with segment length and number of bins, but will be a bit lower if the segment covers a region with less reliable sequencing coverage or mapping (i.e. lower-weight bins).
- The calculation is in the function [cnvlib.fix.apply_weights][3]. It is not one equation; it depends on which data sources are available, mainly whether the reference is paired, pooled, or flat.

Also see `cnvkit.py segmetrics --ci` and `cnvkit.py call --filter ci` for [filtering out potential FP segments][4] by calculating confidence intervals for each segment's mean log2 ratio.

[1]: https://cnvkit.readthedocs.io/en/stable/fileformats.html
[2]: https://cnvkit.readthedocs.io/en/stable/fileformats.html#copy-number-reference-profile-cnn
[3]: https://github.com/etal/cnvkit/blob/master/cnvlib/fix.py#L238
[4]: https://cnvkit.readthedocs.io/en/stable/pipeline.html#filtering-segments
Question: How can I print only 1 column from a BAM file at specific coordinate.

For example:

```
ACGTGACGGACGGGGCCA
TGACGGACGGGGCCA
GGACGGGGCCA
^
```

For example print the column at coord 'x' will return AAA (see example above)

Answer: Maybe something like this?

/package/samtools/samtools mpileup ./your.bam | cut -f 1,2,5

That will give you a chromosome, base position (1-based) and a string similar to what you're looking for, but is case sensitive (strand) and has some other symbols for mismatches, etc.

Check out: http://samtools.sourceforge.net/pileup.shtml
Question: Hello everyone, i wrote this post to know if some of you knew a package on R (or something else) that could create a gene plot such this one :

https://ibb.co/cLm3Tbr

Indeed I wondered if it was possible to `(1)` add the information about the chr or scaffold length and `(2)` order the figure with the longest scaffold on the top and `(3)` add the start and stop coordinates of each scaffold?

It would be very usefull since all scaffold have different length and this information on a plot can be really usefull.

Il already tried gggenes packages but it does not take the scaffold length into account.

Thank you very much for you time and thank you for the package.

Here are the data that could be used for such a plot:

molecule gene start_gene end_gene start_scaff end_scaff strand direction
1 scaffold1 Gene1 64000 68000 1 125000 forward 1
2 scaffold1 Gene2 80000 83000 1 125000 forward 1
3 scaffold2 Gene3 60000 68000 1 80000 forward 1
4 scaffold2 Gene4 20000 28000 1 80000 reverse -1
5 scaffold3 Gene5 22000 29000 1 60000 forward 1
6 scaffold4 Gene6 20000 33000 1 40000 reverse -1
7 scaffold4 Gene7 35000 38000 1 40000 forward 1
8 scaffold5 Gene8 17000 19000 1 20000 reverse -1
9 scaffold6 Gene9 2000 2500 1 5000 forward 1

>dput(tab)
structure(list(molecule = structure(c(1L, 1L, 2L, 2L, 3L, 4L,
4L, 5L, 6L), .Label = c("scaffold1", "scaffold2", "scaffold3",
"scaffold4", "scaffold5", "scaffold6"), class = "factor"), gene = structure(1:9, .Label = c("Gene1",
"Gene2", "Gene3", "Gene4", "Gene5", "Gene6", "Gene7", "Gene8",
"Gene9"), class = "factor"), start_gene = c(64000L, 80000L, 60000L,
20000L, 22000L, 20000L, 35000L, 17000L, 2000L), end_gene = c(68000L,
83000L, 68000L, 28000L, 29000L, 33000L, 38000L, 19000L, 2500L
), start_scaff = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), end_scaff = c(125000L,
125000L, 80000L, 80000L, 60000L, 40000L, 40000L, 20000L, 5000L
), strand = structure(c(1L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 1L), .Label = c("forward",
"reverse"), class = "factor"), direction = c(1L, 1L, 1L, -1L,
1L, -1L, 1L, -1L, 1L)), class = "data.frame", row.names = c(NA,
-9L))

[1]: https://ibb.co/cLm3Tbr
[2]: https://freeimage.host/i/d4imog
[3]: https://ibb.co/cLm3Tbr

Answer: Maybe there is a ready package for this, but we can recreate this using *ggplot2*, here is a start:

library(ggplot2)

ggplot(tab, aes(x = start_scaff, xend = end_scaff,
y = molecule, yend = molecule)) +
geom_segment(size = 3, col = "grey80") +
geom_segment(aes(x = ifelse(direction == 1, start_gene, end_gene),
xend = ifelse(direction == 1, end_gene, start_gene)),
data = tab,
arrow = arrow(length = unit(0.1, "inches")), size = 2) +
geom_text(aes(x = start_gene, y = molecule, label = gene),
data = tab, nudge_y = 0.2) +
scale_y_discrete(limits = rev(levels(tab$molecule))) +
theme_minimal()
Question: I guess it's a simple question.. Could you detail what's the difference between vcf file generated by samtools and from bcftools?

```
samtools mpileup - ref.fa file.bam > file.bcf
bcftools call file.vcf > file2.bcf
```

Answer: Usually these command are used together:

**samtools mpileup** command automatically scans every position supported by an aligned read, computes all the possible genotypes supported by raw reads, and then calculates the probability that each of these genotypes is truly present in your sample.

For example, let's consider the first 1000 bases in Reference Genome file. Suppose the position 35 (in reference G) will have 27 reads with a G base and two reads with a T nucleotide. Total read depth will be 29. In this case, the app concludes with high probability that the sample has a genotype of G, and the T reads are likely due to sequencing errors. In contrast, if the position 400 in reference genome is T, but it is covered by 2 reads with a C base and 66 reads with a G (total read depth equal to 68), it means that the sample more likely will have G genotype.

**bcftools call** command uses the genotype likelihoods generated from samtools mpileup to call genetic variants and outputs the all identified variants.

So, it means, that file.bcf will contain all possible genotypes in the genome, but the bcftools bcf file will contain only sites which were found to be variant.

If you are interested in specific sites that were not called by bcftools, you can break it down into two separate steps.

Do you want to see example vcf files from both commands?
Question: Hi!

I have to undertake a system biology training and I am looking for some online tutorials/courses/guides.

I know a lot of books are available and they deal with it, but does anyone on this forum has notes/instructions/tutorials on the same.

Kindly share them if you have tutorials on the same.

Thank you

Answer: There are a few courses on Systems Biology announced in the Coursera platform:

Dynamical modelling methods for Systems Biology

Introduction to Systems Biology

Network Analysis in Systems Biology

These courses haven't started yet, but usually all the courses in Coursera are good.
Question: How to download all genes inside all gene sets in C6 (oncogenic signature gene sets, 189 gene sets) in Molecular signature database (Msigdb) I want to download the all genes inside 189 gene sets in R in a data frame?

Thanks in advance

Answer: With the [msigdbr](https://cran.r-project.org/web/packages/msigdbr/index.html) R library

```r
library("msigdbr")
library("tidyverse")

df <- bind_rows(msigdbr("Homo sapiens", "C6"))
```
Question: Hello everyone,

I used different tools to convert human gene ensemble ID to ENTREZ but each of them seems to have problems like:
-BioMart does not recognizes all the ensemble ID but barely 30% ,
-BioDBnet ...I couldn't make it working cause I could only select the type of output ,
-DAVID gives me lots of false duplicates.

What do you usually use for this work ?

thank you

Answer: You can use `AnnotationDbi` to convert Ensembl ids. The code snippet converts Ensembl id to Entrez ids.

library("AnnotationDbi")
library("org.Hs.eg.db")
#columns(org.Hs.eg.db) # returns list of available keytypes
data$entrez = mapIds(org.Hs.eg.db,
keys=row.names(data), #Column containing Ensembl gene ids
column="ENTREZID",
keytype="ENSEMBL",
multiVals="first")