cannin/biostars_qa huggingface dataset:
Question: This question was inspired by Peter Cock 's tweet

https://twitter.com/pjacock/status/118012750105546752

Most softwares (C, java...) use a int32 (unsigned or not) to store the length of the chromosome. It isn't enough when the length of a chromosome is greater than INT_MAX or UINT_MAX

# define INT_MAX 2147483647
# define UINT_MAX 4294967295U

So my question is:

Is there any resource where one can find the length of the chromosomes. Something like BioNumbers.

What's the length of the longest chromosome ?

Answer: Second question: *Amoeba dubia*, which if I recall correctly is also known as *Polychaos dubium*, has a genome that is 670 gigabases. While I know nothing about how this is divided among chromosomes, it is reasonable to suspect that a given chromosome would exceed an int32 or even a uint32.

Refs:

1. http://www.sciencedirect.com/science/article/pii/S0169534703003239
2. T. Cavalier-Smith, Introduction: the evolutionary significance of genome size, T. Cavalier-Smith, Editor, The Evolution of Genome Size, John Wiley & Sons (1985), pp. 1–36.
3. C.T. Friz, The biochemical composition of the free-living Amoebae Chaos chaos, Amoeba dubia and Amoeba proteus. Comp. Biochem. Physiol., 26 (1968), pp. 81–90.
Question: Hi all,

In R, I have a data frame with physical position of some genes in bp on the left and intervals (again in bp) on the right. I want to check for each gene, whether it falls in the the interval of the same chromosome.

```
GeneID Chr Position Interval
x 1 18697251 1:3099740-8804450
y 1 19546617 1:3422930-8804450
z 1 3332236 2:2751757-4502486
a 2 3993537 2:6187995-8804450
b 2 9983384 2:4864334-18271005
c 2 11479025 2:8222941-18271005
```

For example: Is gene x on chr 1 comprised in one of the 2 intervals displayed for chr 1? I expect TRUE or FALSE.

I managed to do this giving the gene positions one by one, but the list is very big, I'd like this to be automated.

Any help much appreciated!

Answer: I have a solution with a loop, although it is not very elegant (code not tested, you may need to debug it). I assume that all your data is as factors (this is why I use `as.numeric(as.character())` for instance).

n=length(myData$GeneID)
answer=rep(NA,n) # this will be an array of TRUE/FALSE
for (i in 1:n){
temp = strsplit(as.character(myData$Interval[i]),":")[[1]] # that gives you c("1","309974-880445")
answer[i] = temp[1]==as.character(myData$Chr[i]) #first test chromosome
if (answer[i]){ # no need to go further if the chromosome is wrong
temp = strsplit(temp[2],"-")[[1]] # that gives you c("309974","880445")
pos = as.numeric(as.character(myData$Position[i])) # in case your data is as factor
# Now lets test if position is between intervals.
# I assume here that the interval can be in the "wrong" order and that its ok for you.
answer[i] = (pos < temp[1] & pos > temp2) | (pos > temp[1] & pos < temp2)
}
}

round(sum(answer)/length(answer)*100,2) # give you percentage of rightly positioned genes.
myData$GeneID[!answer] # give you the names of the wrongly positionned genes.
Question: Hello,

How can I export some features (like country and date) from sequences inside a genbank file?

Thank you in advance!

Answer: I could get it using the python script in this [question][1]

[1]: https://www.biostars.org/p/221662/
Question: Hello, I have the following question:
I have an RNAseq dataset (DESeq dataset) that I want to do PCA analysis on. This is how I made the dataset:

dds_norm <- estimateSizeFactors(dds)
dds_norm_vst <- vst(dds_norm)
counts_norm_vst <- assay(dds_norm_vst, normalized = TRUE)

I made a PCA bi-plot, screeplot and pairsplot with the package PCAtools, with the following code:

dds_vst_pca_.1 <- pca(counts_norm_vst, metadata = xp_design_pca, removeVar = .1)
#make scree plot (90% top genes)
scree_all_.1 <- data.frame(percentVar_all_.1)
scree_all_.1[,2] <- c(1:46)
colnames(scree_all_.1) <- c("variance","component_nr")
ggplot(scree_all_.1, mapping=aes(x=component_nr, y=variance)) +
geom_bar(stat="identity") +
ggtitle("screeplot for all timepoints")

#make bi plot (90% top) with PCAtools
biplot(dds_vst_pca_.1,
colby = "time",
shape = "treatment",
pointSize = 5,
legendPosition = 'right')

#find PC's to keep
elbow <- findElbowPoint(dds_vst_pca_.1$variance)

#make pairs plot (90% top)
pairsplot(dds_vst_pca_.1,
components = getComponents(dds_vst_pca_.1, c(1:10)),
colby = "treatment",
colkey = c("e" = "gray",
"m" = "darkgreen",
"n" = "lightgreen"))

Then, I tried to check myself and made another biplot with prcomp() according to the following code:

dds_vst_prc <- prcomp(counts_norm_vst, scale = TRUE)
ntop_.1 <- nrow(counts)*.9
ntop_.1 <- as.integer(ntop_.1)
#do prcomp with lower 10% of variable genes removed
select_.1 <- order(rv, decreasing = TRUE)[seq_len(min(ntop_.1, length(rv)))]
mat_.1 <- t( assay(dds_norm_vst)[select_.1, ] )
dds_norm_vst_prc_.1 <-prcomp(mat_.1)
dds_norm_vst_prc_.1_df <- as.data.frame(dds_norm_vst_prc_.1$x)
percentVar_all_.1 <- dds_norm_vst_prc_.1$sdev^2 / sum(dds_norm_vst_prc_.1$sdev^2)
ggplot(dds_norm_vst_prc_.1_df, aes(PC1, PC2, color = time, shape = treatment)) +
geom_point(size = 4) +
xlab(paste0("PC1: ",percentVar_all_.1[1]*100,"% variance")) +
ylab(paste0("PC2: ",percentVar_all_.1[2]*100,"% variance")) +
coord_fixed() +
ggtitle("PCA plot")

Indeed, these two give exactly the same biplot (same percentages of PC's and same position of the samples).
However, then I looked into scaling, and I see that although for the prcomp() I put scaling to TRUE, for the pcatools method i did not specify, assuming scaling default would be TRUE. However, I found out the default is FALSE! Indeed, when I put scaling to TRUE, I get a different outcome.
How is this possible? What am I overlooking?
Thanks very much for any input!

Answer: `dds_vst_prc <- prcomp(counts_norm_vst, scale = TRUE)`

Is `counts_norm_vst` transposed here? `mat_.1` seems to be, is the inconsistency based on missing transposition?
Question: Hi,

I was wondering if someone could comment on, or point me in the right direction of, considerations when normalising Illumina 450k Methylation data when there are large differences in global methylation status? The experiment I have is where the same cell line is used across 24 samples but there are different treatments and timepoints. One of the treatments is with decitabine, for example, which results in a marked global demethylation seen as a leftward shift in the global beta value profile - see Figure 5A below for an example.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3587326/figure/F5/

I have read around the area a bit, and much of the literature is concerned with exploring differences between cancer/normal or different tissues. This guide from Brent Pedersen was particularly helpful: https://github.com/brentp/450k-analysis-guide

The minfi vignette and associate papers were also useful, but the thing that struck me was the comment in the Dedeurwaerder et al 2014 review also quoted in the minfi Functional Normalisation paper:

> There is to date no between-array normalization method suited to 450K data that can bring enough benefit to counterbalance the strong impairment of data quality they can cause on some data sets

So, am I best off just doing a bare minimum within-array normalisation using, say, the `preprocessRaw` function in minfi and not doing any between-array normalisation at all?

Any comments gratefully received.

Phil Chapman, CRUK Manchester Institute

Answer: Hi Phil,

it seems that your dataset is a perfect example of when to use functional normalization. If there are large differences in global methylation status, functional normalization should be able to keep them while removing unwanted technical variation. Functional normalization is a kosher within-array normalization based only the array control probes. Those, by design, are not associated with the biology of your samples, and therefore global differences in methylation seen between samples should be conserved. In our implementation of functional normalization in minfi, `preprocessFunnorm()` implements as well the 'noob' background correction method (Triche et al., 2013), which improves significantly the downstream analysis results.

Let me know if you have any more questions,

Jean-Philippe Fortin
Question: Hello Biostars Community,

I made a txdb object using:

`mm39.txdb <- makeTxDbFromEnsembl(organism = "Mus musculus")`

and then made the `CompressedGRangesList` :

`txns <- GRangesList(cds(mm39.txdb, columns = c("CDSSTART","CDSEND")))`

----------

**I am trying to figure out how to rename** `CDSSTART` **to** `cdsStart` **and** `CDSEND` **to** `cdsEnd`**.** Help? Please?

The reason for this renaming is so that nothing "breaks" in the package when plugging in the final data files it requires.

I feel like this is such a simple task, but I am not sure exactly to do it... I've probably been trying to figure out this tiny thing for like 2 hours now... I tried exploring `GenomeInfoDb` for something similar to

`seqlevelsStyle(mm39.txdb) <- "UCSC"`

No success, yet...

Thank you in advance!

-Pratik

Answer: You can set the mcols (which is a DataFrame).

mcols(txns, level="within")[, "cdsStart"] <- mcols(txns, level="within")[, "CDSSTART"]
mcols(txns, level="within")[, "cdsEnd"] <- mcols(txns, level="within")[, "CDSEND"]
txns <- txns[, c("cdsStart", "cdsEnd")]

> txns
GRangesList object of length 1:
[[1]]
GRanges object with 528396 ranges and 2 metadata columns:
seqnames ranges strand | cdsStart cdsEnd
|
[1] 1 4878137-4878205 + | 4878137 4878205
[2] 1 4878137-4878205 + | 4878137 4878205
[3] 1 4878137-4878205 + | 4878137 4878205
[4] 1 4878137-4878205 + | 4878137 4878205
[5] 1 4878137-4878205 + | 4878137 4878205
... ... ... ... . ... ...
[528392] JH584304.1 55480-55701 - | 55480 55701
[528393] JH584304.1 56986-57151 - | 56986 57151
[528394] JH584304.1 56986-57151 - | 56986 57151
[528395] JH584304.1 58564-58616 - | 58564 58616
[528396] JH584304.1 58564-58616 - | 58564 58616
-------
seqinfo: 61 sequences (1 circular) from an unspecified genome
Question: Given a set of regions from a CNV analysis, what are the tools folks are using to annotate them, with particular focus on genes, disease associations, and known CNVs? Google was particularly unhelpful, but that may be my limitations with the tool. CNVannotator is prominent, but it is web-only and has a limit of 500 region batches. Any other suggestions?

Answer: We analyze a lot of CNVs both called using tools like Birdsuite or PennCNV or those which are imputed from GWAS using reference panels like that from 1000 Genomes. In terms of tools, I can recommend two steps that we do: the first covers overlap with genes, using Bedtools; the second covers disease associations, using gene2mesh .

From a BED file of CNVs, such as this:

chr22 39378403 39388216 esv2666691 MERGED_DEL_2_106009
chr5 151514804 151518864 esv2666686 MERGED_DEL_2_32905

You take a bed file including all genes in the genome or a subset. You may also use individual exons of genes (contact me for such a file.):

chr1 11873 12227 DDX11L1
chr1 12612 12721 DDX11L1
chr1 13220 14409 DDX11L1
chr1 14361 14829 WASH7P
chr1 14969 15038 WASH7P
chr1 15795 15947 WASH7P
...

Using Bedtools, run the intersection like so:

intersectBed -a cnvs.bed -b refseq_exons.bed -wb

This will give you output like so:

chr22 39387563 39388216 esv2666691 MERGED_DEL_2_106009 chr22 39387563 39394225 APOBEC3B-AS1
chr22 39358280 39388216 esv2666691 MERGED_DEL_2_106009 chr22 39353526 39388783 APOBEC3A_B
chr22 39358280 39359188 esv2666691 MERGED_DEL_2_106009 chr22 39353526 39359188 APOBEC3A
chr22 39378403 39388216 esv2666691 MERGED_DEL_2_106009 chr22 39378403 39388784 APOBEC3B
chr5 151514804 151518864 esv2666686 MERGED_DEL_2_32905 chr5 151338458 151650010 CTB-12O2.1

For your second step on finding links with disease, I have found gene2mesh to be very helpful. It gives links to keywords, but there are some other useful resources including OMIM that may be helpful. In the case of gene2mesh, the following perl script can take the output genes and get the top MeSH terms. In this case, we put the genes from the intersectBed output above into a file called "gene_list.txt". Perl script we modified from the website looks like so:

#!/usr/bin/perl -w
use strict;
use warnings;
use XML::XPath;
use XML::XPath::XMLParser;
use LWP::UserAgent qw($ua get);
my $ua = new LWP::UserAgent;
my $file="gene_list.txt";
open(F,$file);
while(<F>){
( my $gene)=split;
my $getf="http://gene2mesh.ncibi.org/fetch?genesymbol=${gene}&limit=5";
my $response = $ua->get($getf);
my $xp = XML::XPath->new(xml => $response->content);
print "## Top 30 MeSH Terms from Gene2MeSH Associated with GeneID $gene ##\n\n";
foreach my $g2mNode ($xp->find('//Descriptor/Name')->get_nodelist) {
print $g2mNode->string_value . "\n";
}
}
close(F);

And running it on a short list of genes would give:

## Top 30 MeSH Terms from Gene2MeSH Associated with GeneID APOBEC3B ##

Cytidine Deaminase
vif Gene Products, Human Immunodeficiency Virus
Gene Products, vif
HIV-1
HIV Infections
Question: I stated with the wget function:

wget https://github.com/broadinstitute/gatk/releases/download/4.0.2.1/gatk-4.0.2.1.zip

I think my problem deals with the gatk-4.0.2.1.zip file. After I unzip the file I can view the contents; however, no where do I see the **GenomeAnalysisTK.jar** file? I only see the two gatk-package jar files - listed below.

gatk-package-4.0.2.1-spark.jar and gatk-package-4.0.2.1-local.jar

I'm expecting to find the GenomeAnalysisTK.jar file.

Perhaps my instructions are a bit old?

Answer: Now there is a wrapper script called `gatk`, so you just have to call it. Try `gatk` with no parameters, or:

gatk --list Print the list of available tools
gatk Tool --help Print help on a particular tool

The wrapper is a Python script:

#!/usr/bin/env python
#
# Launcher script for GATK tools. Delegates to java -jar, spark-submit, or gcloud as appropriate,
# and sets many important Spark and htsjdk properties before launch.
Question: Hi

What are the best long reads aligners?
I have dna sequencing data generated by minion. Which aligners should I use, I need to detect SNPs and indels downstream

Thanks

Answer: [Minimap2][1] is the currently recommended aligner for long read DNA and cDNA/RNA sequencing. If you are specifically interested in large structural variation you should also try [ngmlr][2].

[1]: https://github.com/lh3/minimap2
[2]: https://github.com/philres/ngmlr
Question: Hi all,

I am comparing two bed files to find correlation between them.now to examine that whether the correlation is occurring by chance or it is meaningful correlation, in order that I have randomized bed file to 1KB upstream and downstream against which I am finding correlation. now I want to find find P value for the same.

A and B are two bed file. A corr B is 0.9 and A corr C is 0.2 . Here how P value will be calculated.

C - randomized B to 1kb

I would really appreciate any help.

Answer: If I understand well, u want to get a pvalue assessing the significance of a correlation.

Using random permutation is a good idea. However you should do the randomization many times in order to have an empirical distribution of the (A,C) correlation. From that distribution, you can then get the significance of your (A,B) correlation.

Good luck !

Edit : The statistical test to use will depend on your distribution. If it is normal (you can use a normality test to make sure of that), then you could compute the mean, standard deviation and pvalue like this (with R) :

#rand_corr = vector of random correlations
#here for testing, 10000 random number with 0.2 mean and 0.2 sd
rand_corr=(rnorm(10000,0.2,0.2))

# pval calculation :

mean <- mean(rand_corr) #mean of random (A,C) correlations
sd <- sd(rand_corr) #sd of random correlation
x <- 0.9 #true (A,B) correlation value
z <- (x-mean)/(sd) #center normalize
2*pnorm(-abs(z)) # return pval

# visual representation :

hist(rand_corr, breaks=30)
abline(v=x, col="red")

However if the normality is not respected (it might not be since you have correlation values that cannot go below 0 or higher than 1), you might need to ressort to other tests.
Question: Previously, I split out a vcf file by chromosome, and for my project, I have combined the X and XY vcf files into a single one. After changing the "XY" chromosome designation to "X" via:

awk '{gsub(/"XY"/, "X"); print;}' Genome_newX.vcf > Genome_newX2.vcf

I'm running into the issue of sorting this new "Genome_newX2.vcf" by position. The idea is that I'll subsequently run the vcf through GenotypeHarmonizer.

Are there any suggestions on how to do this easily? I'm brand new to this style of work, and I'd love some direction on where to read up on it as well. Thank you!

Answer: Edit: 02/23

Just use `bcftools sort` https://samtools.github.io/bcftools/bcftools.html#sort

----------

Original answer with awk-fu:

cat in.vcf | awk '$1 ~ /^#/ {print $0;next} {print $0 | "sort -k1,1 -k2,2n"}' > out_sorted.vcf

It takes a VCF and prints the sorted file including the header.
Question: Hi!

I'm trying to get SomaticSniper to run on Ubuntu 10.04 and I've been following the instruction from the website below:

http://gmt.genome.wustl.edu/somatic-sniper/1.0.2/install.html

However, when I run the command:

sudo apt-get update

I encounter three errors (at the end of the post) along with several successfull hits. Looks like the address http://repo.gsc.wustl.edu is not accessible. I would be greatfull for advise on how to make the installation work.

W: Failed to fetch http://repo.gsc.wustl.edu/ubuntu/dists/lucid-genome-development/Release.gpg Something wicked happened resolving 'repo.gsc.wustl.edu:http' (-5 - No address associated with hostname)

W: Failed to fetch http://repo.gsc.wustl.edu/ubuntu/dists/lucid-genome-development/main/binary-i386/Packages.gz Something wicked happened resolving 'repo.gsc.wustl.edu:http' (-5 - No address associated with hostname)

W: Failed to fetch http://repo.gsc.wustl.edu/ubuntu/dists/lucid-genome-development/main/i18n/Translation-en_US.bz2 Something wicked happened resolving 'repo.gsc.wustl.edu:http' (-5 - No address associated with hostname)

Answer: repo.gsc.wustl.edu is our internal, private APT repo. You should be using apt.genome.wustl.edu. Where did you see repo.gsc.wustl.edu? I'd like to fix the issue.

Edit: It looks like our install.sh is misconfigured to point to our repo.gsc.wustl.edu instead of apt.genome.wustl.edu. We'll work on getting that fixed but in the mean time you can edit /etc/apt/sources.list.d/genome.list to correct the address. After editing do an apt-get update to refresh your APT caches.

Edit: install.sh has been fixed and deployed.
Question: ![**Screenshot: first column (Cell 1) on the X-axis versus all gene names on the Y-Axis.**][1]

I'm trying to create a very simple heatmap from a dataset with no transformation.

Our data table was imported successfully and was defined as follow with 600 genes and 30 cell type (numbered from 1 to 30), the values correspond to gene expression values from microarray data :

```
Gene name Cell1 Cell2
Gene1 120 100.58
Gene2 12.25 1000
Gene3 0.2 10
```

When creating the heatmap, I followed instruction on spotfire website, click on insert > New visualization > Heat Map

Then, on the Y-axis I select the gene names, from here I faced 2 major issues:

1. The heatmap displays only one cell type on the X-Axis, when I click on add, it adds two types of labels, and aggregates with an average value. so 30 cell types turns into 60 columns.
2. I just want to use the values of the dataset, without any aggregations, I don't understand why the aggregation shows up at the bottom of the heatmap and why I can't turn it off? It seems that the option "none" is missing in the drop down menu of the aggregation.

Is there a simple way to just have the gene expression vs cell 1 to 30 on a heatmap in spotfire? (So we should have 600 rows and 30 columns)

[1]: http://i60.tinypic.com/505oo3.png

Answer: If you want a clustered heatmap, you can use the short-wide table in the format you show and the Hierarchical Clustering tool, as mentioned by swbarnes2 above.

Recent Spotfire versions expect a tall-narrow (unpivoted) data table for the standard heat map visualization. If that's what you want, you can unpivot your data outside or inside Spotfire. To do it inside Spotfire,

File > Add Data Table > From (your existing table).

Select Unpivot transformation.

In Unpivot, select Gene column as Passthrough.

Select all Cell columns as Columns to Transform.

Click OK to add the unpivoted tall-narrow table. It will have 3 columns: Gene, Category (Cell), and Value.

Once you have the tall-narrow format data table, add your heat map.

Insert > New Visualization > Heat Map, and visualize the new unpivoted data table.

Set the row selector to Gene; the column selector to Category (Cell).

You'll get a heat map with Genes as rows and Category (Cell type) as the columns.
Question: Hello,

I have some Sanger sequences of HBB gene. I need to find variants related to Sickle Cell Anemia or Beta Thalassemia from these sequences.

Is there any tool for calling variants from these sequences? Please help.

I have .ab1 files.

Thanks

Answer: You can call variants from *.ab1 files with [tracy][1] or the web application [indigo][2], a front end of tracy as described in [this paper][3]

[1]: https://github.com/gear-genomics/tracy
[2]: https://www.gear-genomics.com/indigo/
[3]: https://doi.org/10.1186/s12864-020-6635-8
Question: Hi guys,

I got a metagenome sample and identified several species/families... in there. Now I'm thinking of automatically gathering pathways (just the names though) for each of those species. Preferably this should be done in R. Any Ideas of a Database / API where to download Pathway-Names given the species name?

Thanks,
Phil

to clarify, basically what I need is typing 'Gardnerella Vaginalis' into some KEGG (or any other) API and I retrieve the list of pathways, here just a snapshot

gvg00010 Glycolysis / Gluconeogenesis - Gardnerella vaginalis ATCC 14019
gvg00030 Pentose phosphate pathway - Gardnerella vaginalis ATCC 14019
gvg00040 Pentose and glucuronate interconversions - Gardnerella vaginalis ATCC 14019
gvg00051 Fructose and mannose metabolism - Gardnerella vaginalis ATCC 14019
gvg00052 Galactose metabolism - Gardnerella vaginalis ATCC 14019
gvg00061 Fatty acid biosynthesis - Gardnerella vaginalis ATCC 14019
gvg00071 Fatty acid degradation - Gardnerella vaginalis ATCC 14019
gvg00072 Synthesis and degradation of ketone bodies - Gardnerella vaginalis ATCC 14019
gvg00121 Secondary bile acid biosynthesis - Gardnerella vaginalis ATCC 14019

Answer: https://gist.github.com/RamRS/36e60b08a99aaba36c0fd453bb961e37
Question: I am doing a DGE analysis of a total RNAseq dataset of 2 timepoints (5 reps each). I am particularly interested in looking for changes in expression of 1000 genes.

Currently, I have done the analysis by analysing all genes and then picking the 1000 genes I am interested in.
However, my PI has suggested that I could try doing the the analysis by doing DGE on just the 1000 genes. Theoretically, this should improve the statistical significance since there would be minimal adjustments for multiple hypothesis testing.

Is this an advisable way of doing the analysis? Since differential expression levels are fit to a negative binomial distribution (in the case of DESeq2), wouldn't this just mean most of the 1000 genes I input would end up not being differentially expressed?

Edit: We arrived at the list of 1000 genes as we were interested particularly in genes coding small proteins. Hence, we searched Uniprot for human proteins with a maximum length of 100 amino acids.

Answer: > Is this an advisable way of doing the analysis?

In my opinion, it is not advisable. I would use the entire dataset and then check the p-values of your genes of interest, while being open to other genes that may be statistically significant, too.

Prior to normalisation, you can, of course, rigorously filter your dataset for low-count genes.

Kevin