cannin/biostars_qa huggingface dataset:
Question: Hello,
I use geneSCF since several months without problems to annotate human gene list with sym id against GO_all and KEGG databases.
However, since some days, I encounter problem with KEGG. The problem occur during the step "Retreving gene list for hsa from KEGG" where I obtain a very very short list of genes mapped (3-4 % rather than I obtain before ~70% with the same list).
With the same list with GO_all, there isn't any problem.

Someone have any explanation?
Thanks a lot

Answer: Hi,

Sorry for my late response. Although I follow the tag 'genescf', I did not get any alert from Biostars. Just saw your today's email and updated GeneSCF. Now all the issues are solved. Thanks for your valuable feedback.

**GeneSCF v1.1 (patch release 2) Fixes - 2017/01/13 Release notes**

- Due to changes (January 2017 release) in source databases KEGG and Gene Ontology, GeneSCF database modules has to be updated to work with recent database releases.

**List of solved issues:**

1) Empty output files

2) KEGG not working with GeneSymbols

3) Some organisms from Gene Ontology giving errors or not able to process
Question: Hello guys
I'm working on Rna SEQ data, and I want to get good conceptual knowledge about the analysis of the process.

I used the tool bowtie for indexing my reference genome and then the rest is according to the TUXEDO pipeline.

I don't understand what exactly, "indexing" of the reference genome is and why it had to be done.

I do get the point that there are many genes that can be identified and can be matched to the genome and indexed, But i am not clear as ti why should we do this , how is it done . ( i am new to bioinformatics and with minimal knowledge its so hard to understand the algorithm part of it )

so can anyone put it in simple concepts .

Answer: > Such programs as Maq and Bowtie (Table 1) use a computational strategy known as ‘indexing’ to speed up their mapping algorithms. Like the index at the end of a book, an index of a large DNA sequence allows one to rapidly find shorter sequences embedded within it.

Read more from [here][1]

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2836519/
Question: Hey guys,

I have a list of cDNA sequence ids (e.g., NM_005502.2). How can I get the corresponding Ensembl transcript id, i.e., ENSTXXXXXXXXX.

Thanks

Answer: A really cool ID mapping tool called The Synergizer from the Roth Laboratory might be helpful in this case. Alternatively, you could use BioMart to map your cDNA sequence IDs to ENSEMBL transcript IDs.
Question: My phenotype is not normally distributed. I tried several transformations but none of them seems to improve the normality.
I am now interested in an inverse normal transformation using R.

Is something like the following correct, x being my phenotype?

qx <- qnorm((rank(x)-0.5)/sum(x))

It is based on this paper:
> Yang, Jian, et al. "FTO genotype is associated with phenotypic
> variability of body mass index." Nature 490.7419 (2012): 267.

Histogram before transformation available here: https://imgur.com/a/BwednJB

More information about the raw data:

- min: 10
- max: 750
- median: 54.75
- mean: 86.18217
- variance: 8428.881
- standard deviation: 91.80894

Answer: Okay, yes, here is the page:

Note that they are not transforming the original variable. What they do is the following (for height and weight):

1. build a linear regression model `lm(height ~ age + age^2)`
2. extract residuals from model with `residuals()`
3. transform residuals by inverse norm function `y <‐ qnorm((rank(x,
na.last="keep") ‐ 0.5) / sum(!is.na(x))`

The transformed residuals (squared) are then used in your association test, as follows:

glm(y^2 ~ SNP)

Does that help?

Kevin
Question: Hi all,

What is the difference between GRCh and hg? I have small RNA interaction data with genome like below and they have used hg19 for mapping the reads. I want to use GTF file in Ensemble which is based on GRCh. but I was wondering the hg and GRCh are the same? If not, where can I find the GTF file for the hg19?

chromosome start end interaction ID alignment score of target part strand
chr2 1.9E+08 1.9E+08 human08246 32 +
chr16 5.2E+07 5.2E+07 human04514 32 -
chr3 1.1E+08 1.1E+08 human11115 36 +
chrM	14360	14382	human08122	34	+

Answer: More information to add (courtesy MatthewP): GRCh37 / hg19 / b37 / humanG1Kv37 - Human Reference Discrepancies

In essence: GRCh37 is identical to hg19 on the main contigs (chr1-24), but differ on chrM.

Kevin
Question: Hello,

I am trying to merge multiple VCF files using *bcftools* but it threw an error saying that the file is not compressed.

I want to know if the right command to compress the file would be:

bcftools view -I input.vcf -O z -o output.vcf

Should I index the files too before merging ? If so, which command should I be using ?

Thanks in advance.

Answer: Hello https://www.biostars.org/u/39157/ ,

short answers: Yes and yes :)

Long answer:

Compressing a vcf file can be done in two ways:

1. using `bcftools view` as you show
2. using `bgzip -c input.vcf > output.vcf.gz`

Whether you should index your compressed vcf file or not, depends on what tool you like to run afterwards on that file. Some require it, some not. But the tools that don't need an index file, also doesn't make use of the advantage to have random access to the positions in the vcf file.

There are two indexing format: `tbi`and `csi`. Unfortunately I cannot tell what's the difference between them (maybe other can tell something about it). `tbi`is the standard index format if you use `tabix input.vcf.gz`and `csi`is the standard format if you use `bcftools index input.vcf.gz`. `tabix` and `bcftools`provide both an option to use the other indexing format. At least `bcftools` can work with both index format.

Depending on the size of your vcf file it might be useful to convert them to compressed bcf instead of compressed vcf. bcftools is designed for working with bcf. So in every step it converts vcf to bcf for calculation and back vcf afterwards. This produces a massive overhead. If you have bcf from the beginning it is much faster (In my last test 2-4 times).

fin swimmer
Question: Hi

I'm uncertain about what parameters I'm supposed to feed to FeatureCount.
I have a BAM-file and a GFF file (from Prokka).
The GFF-file looks like this:

```
003__g__Streptococcus_no_23	Prodigal:2.6	CDS	6098	6310	.	-	0	ID=PLMABFGE_03080;inference=ab initio prediction:Prodigal:2.6;locus_tag=PLMABFGE_03080;product=hypothetical protein
003__g__Streptococcus_no_23	Prodigal:2.6	CDS	6340	6504	.	-	0	ID=PLMABFGE_03081;inference=ab initio prediction:Prodigal:2.6;locus_tag=PLMABFGE_03081;product=hypothetical protein
003__g__Streptococcus_no_23	Prodigal:2.6	CDS	6641	7483	.	+	0	ID=PLMABFGE_03082;Name=ykuT;db_xref=COG:COG0668;gene=ykuT;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:O34897;locus_tag=PLMABFGE_03082;product=putative MscS family protein YkuT
003__g__Streptococcus_no_23	Prodigal:2.6	CDS	7514	8083	.	-	0	ID=PLMABFGE_03083;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:P39157;locus_tag=PLMABFGE_03083;note=UPF0340 protein YwlG;product=hypothetical protein
003__g__Streptococcus_no_23	Prodigal:2.6	CDS	8115	8582	.	-	0	ID=PLMABFGE_03084;Name=hmpT_2;db_xref=COG:COG4720;gene=hmpT_2;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:A2RIH3;locus_tag=PLMABFGE_03084;product=Thiamine precursor transporter HmpT
```
The BAM-file is an alignment of my raw FASTQ-files (from a metagenomic sample) mapped against the genes predicted from my assembled bins.
The purpose is to figure out, how many reads are mapped to each gene (the depth/abundance of each gene) and lastly, cluster them together in COGs or some sort of functional categorization.

What should I set my -t and -g values to?

Answer: I found the issue. I was using an alignment of the genes rather than the chromosomes, hence the confusion.
I made a new alignment, and everything looks good now!
Question: Hello,

I am merging 14k vcfs and I get this error:

[W::vcf_parse] INFO 'D' is not defined in the header, assuming Type=String
Error: The INFO field is not defined in the header: D

So, I think I have some vcf that is not well formatted for bcftools, would you know if bcftools has something like vcf-check? I haven't seen it. Then, I could remove this vcf.

Thanks

Answer: here is a command to find the VCF missing the INFO/D in the header

$ find dir -name "*.vcf.gz" | while read F; do gunzip -c "$F" | grep "^##INFO= /dev/null || echo $F ; done
Question: Hi,
I'm trying to use DeSeq2 to calculate gene expression levels.

My count table (all_counts) looks like this

gene cell1 cell2
1 ENSMUST00000000001 712 1613
2 ENSMUST00000000003 0 0
3 ENSMUST00000000010 2 3
4 ENSMUST00000000028 0 0
5 ENSMUST00000000033 0 0

I'm call DESeqDataSetFromMatrix

dds <- DESeqDataSetFromMatrix(countData = all_counts,
colData = data.frame(c('gene', 'cell1', 'cell2')))

But I'm always getting this error

Error in DESeqDataSet(se, design = design, ignoreRank) :
some values in assay are negative

But there are no negative values. I verified this by

> any (!is.integer(all_counts[,2]))
[1] FALSE
> any (!is.integer(all_counts[,3]))
[1] FALSE
> any (all_counts[,2] < 0)
[1] FALSE
> any (all_counts[,3] < 0)
[1] FALSE

Any ideas what I'm doing wrong?

Answer: From the manuals the **countData** must be a numeric *matrix*, from your example **all_counts**, it looks like a *dataframe*.

Try to convert it to matrix:

countDataMatrix <- as.matrix(all_counts[ , -1])

And keep the names of the genes in the rownames (as suggested by @DevonRyan in the comments):

rownames(countDataMatrix) <- all_counts[ , 1]
Question: Hi all,

Just a quick question (as a newcomer to bioinformatics) regarding effect size in differential expression analysis. Why does the field opt for using fold change as a metric of effect size? Fold change doesn't take into account variability, whereas standardized effect size measures like Cohen's d do. So why doesn't the field report effect sizes that take into account variability?

To illustrate an example, say gene X has a mean of 7.20 in condition A and 7.60 in condition B. Fold change for condition B compared to condition A is 7.60/7.20 = 1.05. Say the standard deviation estimates on condition A is 0.09, while in condition B its 0.10. Computing Cohen's d on this, the effect size is somewhere around 4.2, which is a gigantic effect. Fold change and Cohen's d differ dramatically, so why not report effect size estimates that take into account variability rather than fold change?

Thanks,
Mike

Answer: Historical inertia mostly.

Having said that, there's a slowly growing movement within bio-related fields to include confidence intervals in all estimate reports. While those obviously aren't standardized effect sizes, they serve a similar purpose (e.g., you can use them for power calculations, which more people should actually do). Actually, I find confidence intervals more useful than standardized effect sizes, but perhaps that's just me.

I should add that standardized effect sizes are not a panacea. My original background is ion channel biophysics, where we can robustly measure very very very very small changes in an extremely robust manner. I could trivially get Cohen's D values vastly larger than the one you showed, but the results would still be biologically meaningless. That's the thing with relying too much on any single number (this goes doubly for p-values), they're easy to point at and yell, "This is important!", when a finding isn't actually.
Question: For RNA-seq data analysis using DESeq2, a recommended method for batch effect removal is to introduce the batch in the design of the experiment as `design = ~ batch + condition`.

The presence of batch was already known from experiment design and also detected by PCA biplot on the log transformed raw counts. Post DESeq2 how can we assure that the batches have been correctly taken care of? PCA plot on `counts(dds, normalized=TRUE)` did not how any changes.

Correcting the batch using ComBat brings out the desired grouping, however, I understand the output from Combat should not be used with DESeq2. Is there any other suggested tool for differential gene expression analysis post batch correction using ComBat?

Answer: You could include `batch` in the design formula, as you have done, and this will then be 'accounted for' (the statistical inferences will be adjusted for batch) when you perform a differential expression analysis testing for `condition`. `batch` is essentially treated as a covariate in the regression model, just as we do in association studies for `sex`, `smoking status`, `BMI`, etc.

If you then want to directly remove the batch effect for downstream analyses, you would do it on the variance-stabilised or rlog expression values using `limma::removeBatchEffect()`. There was a recent question on Bioconductor:
https://support.bioconductor.org/p/125386/#125387 (follow the link to the DESeq2 vignette)
Question: Hi,

I am trying to run the MuSic analyses on a some cohort data that I have. I am using a 4cpu 2xl EC2 machine. I have been able to apply all analyses up to mutation-relation (MRT). Here the application has been running for several days (running several cohorts in parallel, with sizes ranging between 20 and 200). Non of the runs have completed. Previously, even the most extensive calculations on these cohorts did not take up more than 24h.

Looking at the code I see potentially endless loops such as:

```r
doit=1;
while(doit==1)
{
xi=matrix.sample.v1(x);
if (max(xi)==1) doit=0;
}
```

But there could also be less obvious pitfalls, as R loops - looking for restricted permutations.

Can anyone share from their experience running this? Should I abandon ship?

Thanks,
Ido

Answer: Unfortunately it scales exponentially with the number of genes you're testing. I wouldn't run it on more than 30 genes. There are much better alternatives to this pairwise gene test, assuming you're trying to find pathway-level relations. Try MultiDendrix, MeMO, MUTEX, CoMet, etc.
Question: There are two VCF files that I like to merge them, using GATK or VCFtools. The problem is, they have different chromosomal notation, one has Chr, the other does not. This question could be similar to [this one](/p/13462/)

Is there any quick awk/sed commands that you suggest ?! Also I appreciate if you make comment, which of these two (GATK/VCFtools) is more reliable for this task.

Answer: Give a statistical geneticist an awk line, feed him for a day, teach a statistical geneticist how to awk, feed him for a lifetime... check out this page for find and replace using awk:

http://awk.info/?awk1line
Question: Hi all,

I'm looking to compute the k-th percentile value of some bigWig files (ChIP-seq) across the whole genome. For example, I'd like to know the 99.99th percentile value of a bigWig file across the whole genome.

I'm aware of wiggletools histogram (unfortunately does not provide exact percentile values, which I need) and the BEDOPS answer suggested [in this question][1]. However, the BEDOPS answer only computes percentile values on a per-chromosome basis, whereas I need whole-genome percentile values. I have tried to modify it to compute percentile values on a whole-genome basis to no avail.

If anyone has any ideas, I would appreciate them!

Thanks,

Rachel

[1]: https://www.biostars.org/p/180879/#407367

Answer: Unless I'm missing something about the question, here are a couple possibilities:

Perhaps use a pseudochromosome name for all records:

$ bigWigToBedGraph signal.bw signal.bg
$ awk -vFS="\t" -vOFS="\t" '{ $5 = $4; $4 = $1; $1 = "chrPseudo"; print $0; }' signal.bg | sort-bed - > signal.pseudoChr.bed
$ echo -e "chrPseudo\t0\t248956422" > chrPseudo.extents.bed
$ bedmap --echo --kth 0.9999 chrPseudo.extents.bed signal.pseudoChr.bed > answer.bed

I'm assuming that the ChIP-seq input comes from `hg38`, which is why I create an "extents" file that extends over the longest chromosome in the real `hg38` assembly, i.e., `chr1`. Whatever intervals come out of that `awk` statement, they will never extend past the longest, real chromosome. If you're not using `hg38`, you could change `248956422` to match the longest chromosome in your assembly.

That's a pretty convoluted approach.

Another option might be to just strip all the signal from `signal.bg` into its own text file, and run that file through R.

$ bigWigToBedGraph signal.bw signal.bg
$ cut -f4 signal.bg > signal.txt
$ R
> d <- read.table("signal.txt")
> quantile(d$V1, c(0.9999))
...some answer...

Running `?quantile` in R should give more detail about that command.

Even though I'd love for everyone to use BEDOPS, using R or even Python `numpy.quantile` would probably be faster than rejigging a bigWig into a resorted BED file with a fake chromosome name.

Hope this helps!
Question: I would like to download all the .sra files associated with this experiment (GSE55335). I can do it manually via the website but would like to be able to download them from a script or a terminal. I've downloaded and tried both the sra-toolkit and Entrez direct, but I've not managed to use either successfully for this task. Thank you.

Answer: Have a look at this post
Question: Hi,

One rule in my `snakemake` workflow (in which I run `kraken2`) has a large RAM requirement of ~50 GB. When running the workflow on many samples in parallel (e.g., `snakemake --cores 12 -s workflow.smk`), I get a bunch of errors because the system's RAM is exhausted.

As a result, I've tried limiting the amount of RAM used by `snakemake` by setting a maximum memory usage within the `kraken2` rule:

```
resources:
mem_mb=55000
```

Intuitively, I thought that this would mean that `snakemake` would run each job at a time since running >1 job would exceed 55000 MB (=55 GB) of RAM usage. However, the workflow still fails by running out of RAM despite my system having 64 GB of RAM.

Next, I tried also setting a maximum amount of RAM on the command line using the `--resources` flag:

`snakemake -j12 -s workflow.smk --resources mem_mb=55000`

However, I get a traceback error:

```
Traceback (most recent call last):
File "/home/cfos/.local/bin/snakemake", line 8, in
sys.exit(main())
File "/home/cfos/.local/lib/python3.8/site-packages/snakemake/__init__.py", line 2400, in main
resources = parse_resources(args.resources)
File "/home/cfos/.local/lib/python3.8/site-packages/snakemake/resources.py", line 85, in parse_resources
for res, val in resources_args.items():
AttributeError: 'list' object has no attribute 'items'
```

What am I doing wrong here? Any advice would be appreciated.

Answer: UPDATE: This has been fixed in version [6.5.3](https://github.com/snakemake/snakemake/releases/tag/v6.5.3)
Question: Hello Biostar,

I'm sorry if this is an innapropriate place to ask for it, but I can't find this file by myself and actually don't know if it exist. I'm looking for a file containing the whole systematic classification of RefSeq, I have a table that would look like this in mind :

https://www.encyclopedia.com/reference/encyclopedias-almanacs-transcripts-and-maps/examples-systematic-classification-table

I mostly need the "species name" (latin) associated with the "class". does anyone know if such a file exist and where I could find it ?

Cheers,

Roxane

Answer: You can find this information in the NCBI taxonomy files that can be downloaded from ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump. Please see readme file for more info: ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/new_taxdump/taxdump_readme.txt

`rankedlineage.dmp` contains following columns.

rankedlineage.dmp
-----------------
Select ancestor names for well-established taxonomic ranks (species, genus, family, order, class, phylum, kingdom, superkingdom) file fields:

tax_id -- node id
tax_name -- scientific name of the organism
species -- name of a species (coincide with organism name for species-level nodes)
genus	-- genus name when available
family	-- family name when available
order	-- order name when available
class	-- class name when available
phylum	-- phylum name when available
kingdom	-- kingdom name when available
superkingdom	-- superkingdom (domain) name when available
Question: What books would you recommend for someone with a **bio** and not **maths** background to learn about probability distributions?

I am looking for something that could be easily related to bioinformatics as I work in Human Population Genetics.

Any (free) online courses or good web tutorials are very welcome as well!

Thank you!

Answer: Since you're in the UK, you should check [your favourite bay for second hand things] for the OU courses M248 and subsequently, M343. These will give you a pretty broad intro to probability models and a little bit of statistics (from the former). The latter might get a bit mathematical for you (at the moment). Since you're posting on here, you might be best served by getting [Ewens and Grant][1]

[1]: https://books.google.co.uk/books?id=uZvlBwAAQBAJ&printsec=frontcover&dq=ewens%20and%20grant%20statistical%20methods%20in%20bioinformatics&hl=en&sa=X&ved=0ahUKEwj174-So7_TAhXmCMAKHYJgBrcQ6AEIIzAA#v=onepage&q=ewens%20and%20grant%20statistical%20methods%20in%20bioinformatics&f=false
Question: >1 length=424197 depth=1.05x
ATGCATGCAGTAGCAGATGCAGAGAGACAGATAG
AATAGACAGTAGACGATAGACAGTAGAGATAGAGA
ACGATGATGACCCAGTAGATGACAGTAGACAGATG
>3 length=322465 depth=0.97x
GTACATGGTAGCAGATGCAGAGAGACAGATAGAA
AATAGACAGTAGACTATAGACAGTAGAGATAGAGA
ATATAHATAHATAHAGATTHATHATACAGTATAGAT
>4 length=313463 depth=0.87x
CGTAGATCGTAGCAGATGCAGAGAGACAGATAGA
AATAGACAGTAGACGATAGACAGTAGAGATAGAGA
ATAGCAGTAGCAGTAGCAGATGACAGATGGAGAG
>6 length=285776 depth=0.79x
ATGCATGCAGTAGCAGATGCAGAGAGACAGATAG
AATAGACAGTAGACGATAGACAGTAGAGATAGAGA
TGATGACGATGACGATGGGTAGAACACCAGATGG
>7 length=281883 depth=0.90x
ATGCATGCAGTAGCAGATGCAGAGAGACAGATAG
AATAGACAGTAGACGATAGACAGTAGAGATAGAGA
GTACAGTAGACAGTAGACAGAGGGGGAGATAGGA

I have a fasta contig file which is numbered serially but for some reasons i have to remove some contigs , which resulted in gaps in the contig serial no. as shown above i.e. it became from 1,2,3,4,5,etc to 1,3,4,6,7 and so on.

So,now i want to rearrange the contigs i.e.

1 - 1
3 - 2
4 - 3
6 - 4
7 - 5 and so on till 580 (My total contigs)

AND rename them as

1 - Contig_0001
2 - Contig_0002
3 - Contig_0003 and so on till Contig_0580

I am trying it with awk print f function with increment. But i have not succeeded till now .

Thanks for your valuable time .Can anyone help me in solving this issue.

Answer: Try [seqkit replace](https://bioinf.shenwei.me/seqkit/usage/#replace):

seqkit replace -p '.+' -r 'Contig_{nr}' --nr-width 4 contigs.fasta -o renamed.fasta

Example:

$ echo -ne ">1\nactg\n>3\nACTG\n"
>1
actg
>3
ACTG

$ echo -ne ">1\nactg\n>3\nACTG\n" \
| seqkit replace -p '.+' -r 'Contig_{nr}' --nr-width 4
>Contig_0001
actg
>Contig_0002
ACTG
Question: I'm calling variants using three input/reference files (dedup.bam; genome.fa; and chr17.cds.bed) and creating two vcf files, (1) for gatk.vcf and (2) freebayes.vcf.

GATK returned more variants (POS) and more dbSNPs. Scanning the files you quickly notice the different quality scores between the two files, with far greater range in QUAL (super low and high scores) found in the freebies.vcf file. What factors contribute to the higher variant count in GATK vcf files compared to freebayes?

Answer: Keep in mind that the GATK variant callers are designed to be as sensitive as possible and will therefore include many false positives, so you need to apply some filtering steps after calling to remove those false positives, as described in the GATK Best Practices. It's essentially impossible to answer your question without knowing more about how you did the variant calling in both cases, and what kind of filtering and evaluation you did on the results.

It's also important to understand that QUAL scores are calculated differently by different variant callers, so it's tricky to compare them directly. You'll get more insights from evaluating your results relative to known callsets or truth sets.