cannin/biostars_qa huggingface dataset:
Question: I am trying to do some exploratory bioinformatics on TCGA data using fgsea.

Our lab looks at a specific gene so I was trying to see whether high levels of this gene in TCGA expression data is correlated with enrichment of any genesets. I have been preranking the data using DeSeq2 (and using the F stat as a ranking) and was wondering how I should set up the design.

Because it is a continuous variable I could plug in the scaled normalised counts for this gene straight into the DeSeq2 design or I could split the expression into low/high groups and then run the DeSeq2 to calculate the difference between low/high.

I was wondering whether which of these (if either) is more acceptable?
I assume using the continuous variable makes the most sense but I have only seen it done by splitting the expression into two groups by other bioinformatics.
Is the Wald test with DeSeq2 the most appropriate tool to do this with?

I have run both methods using the hallmark genesets and see very different ranking and similar but slightly different ES results. What are peoples' thoughts?

![NES from HALLMARK][1]
![enter image description here][2]

[1]: https://i.imgur.com/ZY9JNs5.png
[2]: https://i.imgur.com/NO9R1Ii.png

Answer: Personally, I think the low/high stratification isn't ideal because you lose information about the expression of your gene of interest (you're collapsing everything into two values: low or high). I prefer the continuous design (edit: however, please see discussion below; important caveats).

An alternate approach would be to calculate the pair-wise correlation between every gene with respect to your gene of interest (using normalized count values); you can use the correlation coefficients are your ranking. Whether this is "better" than using the deseq2 statistic, I don't know. There are many ways to analyze data and the answer of what is "most acceptable" is not always clear or easy.
Question: Dear All,

I'm trying to download wgEncodeCrgMapabilityAlign100mer.bigWig for human genome 38.

I can get the same file for hg19 from the golden path `ftp://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeMapability/`
I can not find such a file in the hg28 directory
`ftp://hgdownload.soe.ucsc.edu/goldenPath/hg38/encodeDCC/`

What I am missing, any idea?

Thanks!

Answer: http://hgdownload.soe.ucsc.edu/gbdb/hg38/hoffmanMappability/ ?
Question: Hi,

I may have missed something but,

I have 2 ranges objects like this:

suppressPackageStartupMessages(library(GenomicRanges))
#> Warning: package 'S4Vectors' was built under R version 3.5.1

rg1 = GenomicRanges::GRanges(seqnames = c(1,1,2),
IRanges(start = c(1,10,1),
end = c(9,19,20))
)
rg1
#> GRanges object with 3 ranges and 0 metadata columns:
#> seqnames ranges strand
#>
#> [1] 1 1-9 *
#> [2] 1 10-19 *
#> [3] 2 1-20 *
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths

rg2 = GenomicRanges::GRanges(seqnames = c(1),
IRanges(start = c(5),
end = c(15))
)

rg2
#> GRanges object with 1 range and 0 metadata columns:
#> seqnames ranges strand
#>
#> [1] 1 5-15 *
#> -------
#> seqinfo: 1 sequence from an unspecified genome; no seqlengths

Then I want to intesect them by maintaining the adjacent intances in the
original `rg1` into different instances in the result object

The default behaviour merge them in the same instance:

rg_intersect = GenomicRanges::intersect(rg1,rg2)

rg_intersect
#> GRanges object with 1 range and 0 metadata columns:
#> seqnames ranges strand
#>
#> [1] 1 5-15 *
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths

I tried to merge after applying a `findOverlaps` method, because the result
of the FO method only gives 2 hits so I would expect that it also gives 2
instances. However, it also merges the 2 things.

findOverlapPairs(rg1,rg2)
#> Pairs object with 2 pairs and 0 metadata columns:
#> first second
#>
#> [1] 1:1-9 1:5-15
#> [2] 1:10-19 1:5-15

GenomicRanges::intersect(findOverlapPairs(rg1,rg2))
#> GRanges object with 1 range and 0 metadata columns:
#> seqnames ranges strand
#>
#> [1] 1 5-15 *
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths

Created on 2019-03-25 by the [reprex package](http://reprex.tidyverse.org) (v0.2.0).

Answer: Just found that I should be able to do it with the `pintersect`
function from a Pair object. In my example
it works but I am not familiar with this method so I am not sure if it would
be general enough.

suppressPackageStartupMessages(library(GenomicRanges))
#> Warning: package 'S4Vectors' was built under R version 3.5.1

rg1 = GenomicRanges::GRanges(seqnames = c(1,1,2),
IRanges(start = c(1,10,1),
end = c(9,19,20))
)

rg2 = GenomicRanges::GRanges(seqnames = c(1),
IRanges(start = c(5),
end = c(15))
)

GenomicRanges::pintersect(findOverlapPairs(rg1,rg2))
#> GRanges object with 2 ranges and 1 metadata column:
#> seqnames ranges strand | hit
#> |
#> [1] 1 5-9 * | TRUE
#> [2] 1 10-15 * | TRUE
#> -------
#> seqinfo: 2 sequences from an unspecified genome; no seqlengths

Created on 2019-03-25 by the [reprex package](http://reprex.tidyverse.org) (v0.2.0).
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: pass the `-n` flag to [samtools][1]

samtools sort -n inputfile output

http://samtools.sourceforge.net/samtools.shtml

[1]: http://samtools.sourceforge.net/
Question: How to merge unique/non overlapping genes between 2 gene model GFF3 files?

I have 2 gene models:

1. Recent supposedly higher quality gene model
2. Older supposedly lower quality gene model

Some genes though that are described in literature have been removed in the more recent gene model. While they are in the older gene model.
I would like to add these and other "missing" genes back to the more recent gene model, from the older gene model.

But if there is an overlap between genes from the different models, I would like to keep the genes as described in the more recent gene model.

Is there a way to do this?

Answer: `bedtools intersect -v` is one option for identifying the unique genes from the older gene model.
Question: Hello Everyone,

I have a VCF with multiple individuals from multiple populations and I would like to get a summary of the allele frequency spectrum for each population. I know that VCFtools has some nice options for outputting allele frequencies. However my data is from non-model organisms and I think using the reference/derived alleles for calculating allele frequencies is resulting in some serious biases.

So I am looking for two possible solutions:

1. How to calculate the folded allele frequency spectrum (not biased by ancestral/reference allele assumptions) from a VCF file as starting point.
2. How to calculate the allele frequency spectrum using an outgroup species to infer the ancestral allele. Are there any packages out there for this? It does not seem trivial to me to infer which allele is presumed ancestral and to incorporate this into the VCF file to then calculate allele frequency spectrum based on ancestral/derived alleles.

Any comments or ideas are much appreciated. Thank you in advance!

Best
Rubal

Answer: This can easily be accomplished via [VCFlib](https://github.com/ekg/vcflib)

If the end goal is association testing have a look at GPAT: https://github.com/jewmanchue/vcflib/wiki

**EDIT**:

**GPAT++** now supports population summary statistics: https://github.com/jewmanchue/vcflib/wiki/Basic-population-statistics-with-GPAT
Question: I want to store a docker image with the input data, code, software, intermediate results and results. This will be a massive file.

Is there something like the sequence archive for people who want to upload a docker image of their whole study? If not, what is the best/cheapest way to store such a file?

Any related thoughts not directly answering the question are also welcome.

Answer: I get the point of using a Docker image is for reproducibility, but the key selling point of Docker is it's modularity! It seems slightly blasphemous to put *everything* onto a Docker image. Why not just a VM image if you're going to make an offline reproducibility archive of multi-gigabytes?

Instead I would give some serious thought to what genomax said, which is to have a Docker image which automates the process of downloading, decompressing, etc, all the raw public data and turning that into the final result instead. This way, your Docker image would be tiny, you don't have issues with the Docker data and the public data falling out-of-sync if there are ever corrections needed (as a general rule of thumb, there should only be 1 place to download the data from). And of course, you can update your Docker image of 100Mb much easier than updating your Docker image of 100Gb to fix a typo in a script.

As an aside, it seems to be a very popular these days when met with the question of "how can I reproduce this in 10 years time?" people think of the future as being some cataclysmic hellscape where nothing works anymore. Some poor future Bioinformatician slumped over a green and black cathode-ray monitor mumbling about "the wisdom of the ancients" while his buddies peddle bicycles to generate power. All to reproduce the RNA-Seq findings for some 10-year-old study.

As someone typing this while playing Pokemon Red via Gameboy emulator on his phone, i'd say the chances of bad code still working in 10 years from now is fairly high, so long as the code was fairly popular at the time :)
Question: Hi,

I have got my first Nanopore sequencing data and the first step was to see if the data is good. Has anyone has any experience with this kind of data and can tell me how to interpret the results.

The whole report can be downloaded [here][1] (not sure how to post it here). Allin all it looks quite good to me, but what I'm not sure about are the two images attached here. These are the per base QC and sequencing content. It seems that the beginning of the reads is not good, but this can be probably trimmed by removing adapters etc. But Is it really that i need to remove the first 1000 positions? This seems a bit extreme.

About the sequencing content, I don't know where to begin. This looks consistent with the first image, where the quality is not good, and if removing these positions, it should get better. but is it ok for the two pairs TG and AC to be so much apart like that?

EDIT:

I have done both pycoQC as well as the minionQC run (R). They show similar results. It can be downloaded from [here][2]
a few of the images are also attached.

![Basecalled reads PHRED quality][3]
![Output over experiment time][4]
![Basecalled reads length vs reads PHRED quality][6]
![Channel activity over time][5]

To me it looks as if the run was past its prime after ~50h. There is no real gain of new reads afterwards. As I'm usually work with mRNA-Seq, I'm not sure how to call the PHRED Quality. I can see that the most reads are of Q ~=8 and most of them are short (which is also expected). But all in all can this run be classified as good?

thanks

![per base QC ][7]
![per base Seq content][8]

[1]: https://datashare.biochem.mpg.de/s/CL7lvs0b1wMSlZj
[2]:https://datashare.biochem.mpg.de/s/lkhPxr417F4eI6n
[3]: /media/images/f4b05e56-2eab-41f0-8b11-46dc209b
[4]: /media/images/cfd343fd-d3ad-41e6-8d8a-6aa6c81c
[5]: /media/images/5fb9d993-e5b2-411e-94b6-40df8be8
[6]: /media/images/f880adb8-691a-4a7a-bfdb-7a185cef
[7]: /media/images/0f7bf463-05c5-467c-8b94-24c6abd4
[8]: /media/images/84f7f372-dfb1-438a-a42f-0af1e6b0

Answer: a few of your reads are very long and those skew and alter the plots.

there is also no binning for the first 10 bases then it is binning into huge bins subsequently, again makes the plot misleading

do not remove the beginning of the reads for QC reasons, that is rarely an advisable course of action.

I would filter the reads to a more manageable length, say 15K and rerun the QC analysis
Question: I am using esearch query as `$query = "SS1G_03709+AND+gene[filter]";`, but it gives me all (gene+ mRNA+genome sequences). What filter do I need to use so I only get gene sequences in my search? I tried a few filters from [here][1], but couldn't find anything to limit my search for genes.

[1]: https://www.ncbi.nlm.nih.gov/books/NBK49540/

Answer: One way (sequence truncated for brevity):

$ esearch -db nuccore -query "SS1G_01676 [GENE]" | efetch -format fasta | grep ">" | grep -v "genome" | awk '{print $1}' | epost -db nuccore | efetch -format fasta
>XM_001597432.1 Sclerotinia sclerotiorum 1980 UF-70 hypothetical protein partial mRNA
ATGGCGCCCAAATTTTCGGAAGACGAAATTGATGATTTGATATATTTTGCTCGTATTGGGGACAATGATG
AATTCGAGAAATTGAGGGAGGAGTTATGTAAAAGGGAAGGATGTTCGATTGCTGAACTGCTAGAGACTGC

$ esearch -db nuccore -query "SS1G_03709 [GENE]" | efetch -format fasta | grep ">" | grep -v "genome" | awk '{print $1}' | epost -db nuccore | efetch -format fasta
>XM_001595570.1 Sclerotinia sclerotiorum 1980 UF-70 hypothetical protein partial mRNA
ATGCATTTCTCAACTGCAAAAACGCTTCTTCCTCTCGCAGTTCTAGTTTCCTATACCACCGCTCAAACAA
CAGCTGCAGCACCACCTGTTGCTAGTGCTCCTACAGGCGGCACTTCTAGTACTTGTCTCGGACAAAATGT

@vkkodali has much neater ways of finding this info in that answer. Taking some inspiration from one of the command there you could ust do

$ esearch -db nuccore -query "Sclerotinia sclerotiorum 1980 [TITLE]" | efilter -molecule mrna | efetch -format fasta > s_sclerot.fa

to get all of them at one time.
Question: Hi, everyone:

I found that there was no golden standard method for single cell RNA-seq subgroup. There were many papers showed different methods for this analysis. Now I wanna try these different clustering methods to my data. How to evaluate them using R, python, etc. I don't want to do biological experiments~

Generally, is there a golden standard (methods, books, etc.) for evaluating clustering methods by computer in machine learning area?

Thanks very much !

Answer: There is no standard method for evaluating clustering methods other than comparing to some ground truth when it is available. There are however several measures of clustering "quality" e.g. the silhouette (in R silhouette). See for example the python sklearn page on evaluating clustering.
Question: Hey I have a fasta file which contains more than 80 sequences. I want to take each sequences in different files and the file name should be header of corresponding sequence.

Answer: curl https://raw.githubusercontent.com/gouthamatla/fasta_File_Manipulation/master/SplitFastaFile.py | python -
Question: Hey guys!

Just a newbie in bioinformatics. I have a problem about integrative bioinformatics: how to combine those database to make my queries more efficiency?

Is there any browser or interface which can help me to search multiple databases while only input one query?

Answer: As Sean said, there is no such magic box for all queries for all databases, nevertheless you could have a look at bio2rdf to get an idea of what is a Integrative database.

http://www.ncbi.nlm.nih.gov/pubmed/18472304

J Biomed Inform. 2008 Oct;41(5):706-16. doi:
10.1016/j.jbi.2008.03.004.

Bio2RDF: towards a mashup to build bioinformatics knowledge systems.

Belleau F, Nolin MA, Tourigny N, Rigault P, Morissette J.

or look at the tools provided by the NCBI http://www.ncbi.nlm.nih.gov/books/NBK25500/
Question: Hello everyone,

I'm comparing 3 sets of variants from 3 different NGS pipeline runs, where all 3 operate on a common set of 150 samples. I see a lot of difference in the variants found in each run.

The three runs are structured as follows:

1. The 150 samples are run as a single batch
2. The 150 samples are run as part of a larger cohort (total cohort size = 250 samples, say)
3. The 150 samples are split into 5 batches of 30 samples each

For #2, GATK SelectVariants is used to extract variants found in the 150 samples. For #3, CombineVariants is used to combine the VCF files.

When I look at a venn diagram of these 3 sets, I see only around a 40% overlap in variants. For reference, 100% is the set of all unique variants discovered across all 3 runs.

To exclude pipeline quirks (AKA "this is the default behavior"), I compared 2 runs that were run 6 months apart on the exact same 25+ samples, and the variants discovered were identical. So, we can confidently say that only the cohort size difference could have caused this gap in variant discovery.

Can we discuss what could be the case here please? I wish to understand why I see 3/5 of the dataset not being called in at least one of the runs.

EDIT: These 3 runs are only computational NGS pipeline runs, they're not sequencing runs. In other words, I'm using the same BAM files across the board.

Answer: Correct me if I misunderstood your question, but it seems like you are splitting or combining samples into various sized groups and then running those through GATK.

From the supplemental material of GATK on haplotype calling:
> Whereas the initial phase of the algorithm is run per sample, the second stage combines the genotype
likelihoods over all samples in order to determine the most likely alternate allele frequency in the cohort.
The likelihood for a given set of genotype assignments at a given frequency is simply the product of the
genotype likelihoods for each sample given that sample’s assigned genotype (Box 1). We then apply a
population genetic prior to the allele frequency likelihoods based on θ, the population specific
heterozygosity, and choose the most likely allele frequency and associated genotype assignments. The
variant quality score for a polymorphic call is given as -10*log10(probability that the site is actually
monomorphic).

Most variant callers use some kind of Bayesian modeling for within-group variation, since being able to model within-group variation significantly increases the likelihood that the remaining variation is due to inter-group variation (usually treatment vs control, tumor vs normal in most experimental designs). Splitting the samples into artificial cohorts of sizes that are not true to the actual data will affect the results of this calculation substantially.
Question: Hi Guys,
I am new to the forum to come directly with a question / problem.

I want to write a program that counts the patterns like this.
I have given such sequences as an example:

> - A: AAAAAAAAA
- B: AAAAACAAG
- C: GAAAACGAA
- D: AAAAAAAAT
>

A patten is a group of four letters of the four groups. As an example, the first Patten would be "AAGA". Now I wonder how many times that happens in my entire Sequenze. So far I have only found programs that are looking for a given pattern. Does anyone have an idea how to program this or do you already know where this was done?

I'm very thankful for your help!

Best regards

Answer: What you are looking for is 'kmer frequency'. There are lots of pre-existing tools to calculate this.

Our very own Alex Reynolds has this github repository for instance:

https://github.com/alexpreynolds/kmer-counter

Kmer frequency will give you the number of occurrences of ALL strings of `k` letters long. If you're only interested in particular kmer occurences, you can filter the data you get back. You may even be able to specify which kmers you want counted specifically.

The above is not a python solution, but there are plenty out there that are which you should be able to find easily now you have that key terminology.

[Here's][1] a tutorial you could follow if you wanted to code something yourself though.

[1]: http://claresloggett.github.io/python_workshops/kmer_counting.html
Question: Dear all,

referring to the batch correction methods for scRNA-seq, would you have any preference and/or comments ? among possible choices :

-- MNNCorrect, as outlined in SimpleSingleCell workflows :

https://bioconductor.org/packages/release/workflows/html/simpleSingleCell.html

-- ZINB-WAVE :

https://bioconductor.org/packages/release/bioc/html/zinbwave.html

-- HARMONY :

https://www.biorxiv.org/content/10.1101/461954v2

-- SCTransform :

https://satijalab.org/seurat/v3.0/integration.html

thanks a lot,

bogdan

Answer: From experience, SCTransform does not perform well unless the majority of the cells are of the same type. It will force true unique populations together with a heavy hand, whereas MNN is much more orthogonal in its changes. Seurat even has a wrapper around [fastMNN](https://htmlpreview.github.io/?https://github.com/satijalab/seurat.wrappers/blob/master/docs/fast_mnn.html).

Haven't tried the other options though, so can't speak to them.
Question: Hi
I have a data frame with multiple columns indicating SNPs ID, chromosome number and position of the SNP (BP):

df1 SNP CHR BP

2	2	rs9391724	6	31320795	5.41323	2	3.83E-07	6.417103
3	3	rs147949474	4	100738958	5.38602	3	5.74E-07	6.241012
4	4	rs3819285	6	31322742	5.32917	4	7.65E-07	6.116073
5	5	rs116548975	6	31320833	5.29763	5	9.57E-07	6.019163
6	6	rs138441284	4	100338276	5.26365	6	1.15E-06	5.939982
7	7	rs9391848	6	31320145	5.24588	7	1.34E-06	5.873035
8	8	rs181964803	6	31319900	5.1553	8	1.53E-06	5.815043
9	9	rs1377457	11	61144652	5.06155	9	1.72E-06	5.763891
10	10	rs9266070	6	31319618	5.03007	10	1.91E-06	5.718133
11	11	rs28835675	19	6807762	4.99524	11	2.11E-06	5.676741
12	12	rs9266066	6	31319525	4.98046	12	2.30E-06	5.638952
13	13	rs79215153	4	100398008	4.97751	13	2.49E-06	5.60419
14	14	rs368752130	6	31320017	4.90517	14	2.68E-06	5.572005
15	15	rs28873729	19	6807814	4.86929	15	2.87E-06	5.542042
16	16	rs5771860	22	49203073	4.8598	16	3.06E-06	5.514013
17	17	6:32451822	6	32451822	4.85825	17	3.25E-06	5.487684
18	18	rs12791961	11	61152028	4.80564	18	3.44E-06	5.462861
19	19	rs35957957	6	31320064	4.77551	19	3.64E-06	5.43938
20	20	rs10897155	11	61141164	4.77384	20	3.83E-06	5.417103
........

Then I have another data frame with various gene names with indicated the start and end of that particular gene:

df2 gene CHR. start end

2	ADCY9	16	4003388	4166186
3	ADORA2B	17	15848231	15879060
4	ATP2B4	1	203595689	203713209
5	C6	5	41142336	41261540
6	CD36	7	79998891	80308593
7	CD40LG	X	135730352	135742549
8	CDH13	16	82660408	83830204
9	CFTR	7	117105838	117356025
10	CR1	1	207669492	207813992
.......

I want to subset the first data frame and get only the SNPs that fall into the gene coordinates based on the rows values from the second data frame (that has to be done for each gene)
I wrote something like this but the process is taking long time so I suppose there is something wrong:

for(i in 1:nrow(df1))
for(i in 1:nrow(df2)){
subset(df1, df1$BP > df2$start_position & df1$BP < df2$end_position)
}

Any help highly appreciated, thanks

Answer: We can use ***data.table::foverlaps***, see example:

# example input
snp <- read.table(text = "snp	chr	bp
snp1	10	111
snp2	11	222
snp3	12	333
snp4	Y	444
", header = TRUE, stringsAsFactors = FALSE)

gene <- read.table(text = "gene	chr	start	end
gene1	1	10	200
gene2	11	100	300
gene3	12	300	350
gene4	X	1	100
", header = TRUE, stringsAsFactors = FALSE)

library(data.table)

# both genes and snps must have the same keys. Adding bp as start and end columns.
snp$start <- snp$end <- snp$bp

# now convert to data.table object with the same keys.
setDT(snp, key = c("chr", "start", "end"))
setDT(gene, key = c("chr", "start", "end"))

# and merge on overlap
foverlaps(snp, gene, type = "within", nomatch = 0L)
# chr gene start end snp bp i.end i.start
# 1: 11 gene2 100 300 snp2 222 222 222
# 2: 12 gene3 300 350 snp3 333 333 333