text:
STAR-long parameters for aligning RNA ONT reads to genome
Are there any suggested parameters to align ONT reads to the reference genome using STAR-long? For now, I used the parameters suggested here, but I noticed a weird behaviour.
I have RNA reads (D. melanogaster) from R7 and R9 flowcells, separately. I only selected to analyze 2D reads in pass category.
I have, respectively, 113249 reads for R7 and 40318 reads for R9. I aligned those reads and get (only!) 150 uniquely mapped reads for R7 data and 8017 uniquely mapped reads for R9 data. I tried to run again the same command on a different server with fresh compilation, but the output file is consistent with these 150 reads.
However, if I align the same with GMAP, I get 78016 uniquely mapped reads for R7 and 33523 uniquely mapped reads for R9, so I suspect that something went wrong in the alignment run.
I am aware that the two mappers behave very differently, STAR-long being  more precise and preferring to report mappings of fewer reads but at better loci, and GMAP being overall less precise, trying to map the most of the reads but at not-so-good loci.
I was wondering if some of you had experience with this and could suggest me the best parameters for RNA reads from ONT?

I guess my reply would be more appropriate as a comment - but I can't. I'm not sure if that setting of SE makes sense, but whatever. This is not necessarily a solution for the behaviour you see, but I would suggest trying https://github.com/lh3/minimap2 for cDNA-seq alignment.

At my hand, STAR-long doesn't work well with noisy reads. It worked for iso-seq mostly because the error rate of iso-seq reads is much lower.

I've had great results using minimap2, particularly when combined with a pre-treatment of Canu for error correction (using minimap2 for the read-to-read mapping):
# correct reads
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
   genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
   -p 4T1_BC06 -d 4T1_BC06 \
   -nanopore-raw workspace/pass/barcode06/fastq_runid_*.fastq
# align reads
~/install/minimap2/minimap2 -p 10 \
   -a ~/db/fasta/mmus/ucsc/mmus_ucsc_all_cdna.idx \
   -x splice <(pv all/4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  samtools sort > 4T1_BC06_all_vs_mmusAll.bam

Update: the most recent nanopore basecaller, combined with the most recent minimap2, seems to now do a decent job with mapping, and pre-correction of reads no longer seems necessary. More recently I've been using LAST to map to the transcriptome, and minimap2 to map to the genome. Minimap2 has the ability to use a homopolymer-compressed genome index, which means that the most common consistent error for nanopore reads (i.e. misjudgement of homopolymer length) will not influence the mapping rate.

Yep I also ended up using minimap2, thanks. At the time of this question, I was just trying to check if I was doing some nonsense with parameters for STAR, or if its behavior was simply due to its algorithm. For reads correction, aren't you afraid of altering the indel composition of your reads if you use Canu?

No... why would that be a problem? Canu is correcting errors based on the consensus of overlapping reads; I don't see why/how that correction would cause more problems than it solves.


text:
What is the best distribution to model the FPKM values from normalized RNA-Seq data?
I know that the discrete raw counts from the RNA-Seq data are usually modeled by a negative binomial or a Poisson distribution, but what I am working on are the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values which represent a transformation of the raw values. They are continuous except many actual 0s. What is the probability distribution FPKM values are likely to follow, or are assumed to follow? Even after taking logarithm twice (after replacing 0s with 1, of course), the data does not follow a normal distribution since there are a lot of 0s.

I think you may want to look into zero-inflated distributions?

Could you show or provide some sample data to explain with? What have you read about limma, DESeq2, DESeq ?

FPKM/TPM values are generally log-normally distributed. Reference : Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels

Sorry, I added a reference. I'm not sure why mRNA levels follow that distribution but it has been established. As they mention in the paper, it's a fairly common distribution in nature. On the other hand, 3' capture methods of single-cell RNA seq such as DropSeq usually yield a negative binomial distribution, presumably because of the high drop out rates of this method of transcript capture.

Well done. The answer looks WAY better with the reference. If you would like to push it even further. This type of information you have written in the comment you can add to answer as well to give a complete picture.

If some data follows a log-normal distribution, does it mean that it would follow normal distribution after you take the log of the values? That's not what I am seeing...

It should...what does expression look like for a gene you know is expressed in your sample after you log normalize your data?

Cuffdiff is one the tools that deals FPKM data for analysis of differential expression.By default they use the negative binomial with some condition specific parameter to avoid biases. There are some other models also available in Cuffdiff. Check the following documents for more details,
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html#library-types
https://www.reddit.com/r/bioinformatics/comments/3bx3em/fpkm_vs_raw_read_count_for_differential/

FPKM data is not composed of integers, it is interesting that they are using negative binomial distribution.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Counting a specific consecutive character with its occurrence position and length
I have a sequence file and want to count consecutive character "N", with the tandem's position of occurrence and its length. Say I have a file named mySequence.fasta like this:
>sequence-1
ATCGCTAGCATNNNNNNNNNNNNNNCTAGCATCATGCNNNNNNATACGCATCACANNNNNNNNNCgcatATCAC

and anticipated output should be like this:
Position 12 N 15
Position 38 N 6
Position 56 N 9

Kindly help me to solve this by awk or sed providing my file name mySequence.fasta

Why are you tied to awk and sed?  There are many tools in programming languages like python, perl and java that are designed for this.

I couldn't come up with a one liner but have a two liner instead...
To get the start positions of the N repeat:
grep N mySeq.fasta | grep -b -o N | awk -F: '{print $1}' |awk '$1!=p+1{print $1}{p=$1}'
12
38
56

To get the length of the repeats:
grep N mySeq.fasta | awk -F '[^N]+' '{for (i=1; i<=NF; i++) if ($i != "") print length($i)}'
14
6
9

With thanks to stack overflow.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
checks for spike-in sequence controls
I would like to know what do people verify when designing/using spike-in controls, to be used in sequencing experiments (mainly Illumina). So far I came up with this list:

Does it align only to a given genome synthetic reference?
Does it contain G-quadruplexes? E.g. QGRS Mapper
Does it contain known Illumina error motifs? (Do this matter anymore for recent platforms/chemistries?) See this article
Does it contain other motifs known to be polymerase inhibitors? Ref

Anything else?

Disclaimer: I'm a developer for http://www.sequin.xyz.
Sequin is a new set of spiked-in controls for next-generation sequencing, and that includes Illumina. We design controls for RNA-Seq, genome sequencing, metagenomics etc.
Please study our papers if you want more details.

Reference Standards For Next-Generation Sequencing by Hardwick should be a good start.

Here's my list:

Measured counts against input concentration. Pearson's correlation, spearman correlation and regression slope.
PCA analysis
Whether the synthetic aligned to the synthetic reference
Spiked-in amount (e.g. dilution)
Diagnostic statistics
ROC Curve
Detection limit

(1) is important because the controls should give you correlation = 1, and slope = 1 in a perfect flawless experiment. This is a sample output (you can read more about it in our papers):

(2) This is our sample PCA (our paper has the details)

(3) is also important but you can only do that in simulations.
(5)

(6)

(7)

Does Anaquin/Sequin allow for the design of new spike-in controls that don't align to a set of references?

@719016 Sequins assume the spiked-in reads aligned to a synthetic reference.

Does it allow for the design of new spikeins wrt a new synthetic reference?

Does it check if the set of synthetic reference is unique against other Genomes and unique against itself?


text:
Run cuffcompare in strand-agnostic mode
Is there a way to run Cufflinks' cuffcompare in a strand-agnostic mode?
I would like to do this because I have some RNA-seq datasets derived from an unstranded run, that should be compared to a reference transcriptome. I think that it makes no reason to have some transcripts in the experimental dataset classified with transfrag code tag x and s as compared to the reference if the input dataset was not derived from stranded protocol.
From manual:

x : Exonic overlap with reference on the opposite strand
s : An intron of the transfrag overlaps a reference intron on the opposite strand (likely due to read mapping errors)

Conceptually, does it make sense?
If so, and if the run by cuffcompare cannot be changed to a strand-agnostic manner, how should I interprete those transcripts classified as x? I was thinking of merging them to o and/or c .
Again, from manual:

c :   Contained
o :   Generic exonic overlap with a reference transcript

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How can I work out gram staining information from the species name?
I have a very long list of prokaryote organism names. I would like to know their gram staining (gram positive or gram negative) information.
Please let me know if there is any database that keeps this information.
I am also trying to find a program that can find gram staining information via organism name only as my organism list is long.
Edit: for each species the NCBI taxonomy lineage is also available

You only have their names, or you also have their genomes?

I have only their names, I can able find their taxonomy lineage using NCBI taxonomy database if that helps me to find the gram staining.

The GOLD database can be used to retrieve a list of species names according to their Gram staining: go to the "Advanced search" page and select the either "gram-" or "gram+" under "Organism Fields" -> "Gram Stain". Then on the results page click on the number under "Organisms". On the results page you can select to have 1000 results per page, meaning that you would need to do some manual work.
I should add that as far as I know this information is based on what people submit when they register their genome projects, so I wouldn't expect this list to be curated or 100% accurate.
If you have access to the genome of the species you are interested in (you could download a representative strain for each species using ncbi-genome-download), you could use traitar, which can predict the gram staining, among other interesting broad phenotypes.

Thanks for your help. I was thinking about downloading GOLD database data and use a program to find my required information. But I am unable to download the data. So can you please help me to do that.

Is there any way to query the GOLD database programatically?

No, their usage policy is in fact quite strict: "Downloading GOLD metadata by any means including web searches, downloadable reports, screen scraping, etc., is strictly prohibited without prior permission, if the purpose is to dissect, alter and store the GOLD data in a database (relational or text) and distribute to third parties/general public/collaborators."

AMR R package can assign gram staining based on taxonomy. See mo_gramstain(x, language = get_AMR_locale(), ...).

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Ultimate reproducibility in R?
I'm looking for a convenient and reliable way to make an R analysis reproducible, either at different times or across collaborators.
Listing the package versions or a sessionInfo() output is not very helpful, as it leaves the work of re-creating the environment entirely to you.
As I see, in order to make analysis reproducible, the following has to be the same:

the R environment (R, repos (Bioconductor), packages, and maybe RStudio)
data
code
seed for random processes: set.seed()

I'll post an answer, but would like to hear about other options.

tl;dr: There is currently no way to accomplish this in R. packrat comes close but ultimately falls short, because it only supports loading one version per package, but different packages might share the same dependency, but require different versions. R needs something like Ruby’s Bundler. I’m planning to work on it, but it won’t work with packages due to fundamental architectural limitations.

You're right about trying to get the R environments to be as similar as possible, as easily as possible. It's not always necessary to have the exact environment to get the same results, but it definitely helps if you find yourself struggling to replicate results.
General tips for ensuring reproducibility, not specific to R
See Sandve et al., PLoS Comp Bio, 2013. They mention a number of things that you do, like recording seed values and session info. I'll highlight a few points:

Environment: Include a small script that lists all the necessary dependencies

This should include all version numbers of languages and packages used
Better implementation: make it a script that automatically installs dependencies

Data: Share raw and/or processed data where possible

At the very least, details on how to obtain the data you used
This can be an accession number in a database, or a link to the source material
Better implementation: again, make it a script that automatically downloads the raw data you used

Code: Comment and document your code

People (yourself included) need to understand what you've done and why
Additionally, your work may be reproducible, but it may not be correct. If you (or anyone else) attempts to fix a problem in your code, they will have a difficult time doing so if they can't understand what your code is doing

Code: Use version control

Services like GitHub and BitBucket make it simple to share your exact work with collaborators
Version control also ensures that you have a history of all the changes done to your work

No manual manipulation of data

If something is done to your data, write it down

R-specific tips

The packrat package, as you mentioned in your answer
Using a number of custom functions for analyzing your data? Compile those functions into an R package

Allows others to easily load/unload the exact functions you used through tools like CRAN, BioConductor, or devtools::install_github()
Makes it easy to see whether these functions are working as intended, since they are separated from other parts of the analysis

Make unit tests for your work with the testthat package

Helps ensure that your code is working as intended, and gives others confidence in your work as well
Write the tests using testthat::test_that(), and check that they pass using R CMD check or devtools::check()
This works especially well with R packages
This also works especially well with TravisCI, since you can automatically have your code get tested after you change it and make a commit

Consider Docker images for your environment

This bundles everything about the R environment together into one easily executable image (and it's not just limited to R)
These images can be easily shared through Docker Hub
Note this might not always be the best solution, since few people know Docker and it's not always simple to use. Docker also requires admin privileges, which may not be possible for work on computer clusters
There is lots of current discussion about Docker in science. See Boettiger, 2014, for an example, using R

If you use Docker, Bioconductor has images too.

Another option I'm trying at the moment is the combination of snakemake, conda.io and knitr.

Conda is a "Package, dependency and environment management for any language" and could be used to set up the R-environment as well as other software. Note, most bioinformatics software can be found in the related bioconda project.
Snakemake is a  workflow management system, which can be used to automate the preprocessing of your data in a reproducible manner (also it works nicely with conda since a conda environments can be specified in the Snake-Makefile)
Knitr for dynamic report generation with R. Well suited to document and reproduce your R workflow into reports, ticking your code and seed boxes.

how is snakemake working nicely with conda, as respect to more "traditional" systems like Makefiles?

@mgalardini I read it as "is installable and works if installed via conda".

There are dedicated directives for conda in the Snakefile. I changed the text to make it more clear. Also have a look here for further info

checkpoint is a system for ensuring you get packages of the same time / version. However, I confess that I haven't used it and thus can't vouch for it.
My feeling is that the rest will have to be handled outside this: code and data can be versioned in git (and dat). A whole solution might have use something like docker as others have suggested.

One way to achieve this is via packrat, a dependency management system developed by RStudio.
Setup:
1) Create project
2) Create script in project directory to start RStudio with specific R version:
export RSTUDIO_WHICH_R=/path/to/R
rstudio

3) Bring project under packrat control. Install bioconductor, install packages...

Every time you work on the project, run the above script in project directory.
https://github.com/rstudio/packrat/issues/342#issuecomment-315059897


text:
What will be an appropriate mathematical distribution for SNP data?
I found that several papers describe SNPs as a binomial distribution with the probability of "success" equals to minor allele frequency.
However, in my experiments, when I generate SNP array following this distribution, the simulation results behave very differently from the SNP array generated following the random mating procedures. I feel like there is some missing information about the description of a binomial distribution.
I wonder what will be a better mathematical description of the SNP array data? Besides a direct answer, any thoughts, or any suggestion of papers are also highly welcome.
EDIT:
Thanks for these suggestions, after reading relevant materials, I think what I observe is about the Hardy-Weinberg equilibrium affected by evolutionary influences like mutation and migration. I will continue to look and list findings here if I find anything interesting. Please feel free to suggest some mathematical descriptions.
Also thanks for the suggestion of tools. They look very helpful, but I really hope to understand some mathematics behind it.

Are you familiar with linkage disequilibrium and Hardy-Weinberg equilibrium? Looking up review papers on those might give more insights.

Re the edit: given the very low mutation rates in eukaryotes, mutation effect on HWE in SNPs should be negligible - or at least, much lower than the effect of genotyping errors.

Simulating genotypes with realistic correlation structures is indeed not so simple, and there's quite a few papers dedicated entirely to that (e.g. https://bmcgenet.biomedcentral.com/articles/10.1186/s12863-015-0173-4). Also, DEPICT (https://data.broadinstitute.org/mpg/depict/index.html) comes with a number of simulated GWASs to generate the nulls, so that's a simple way to obtain some ready-made fake data.
On the other hand, indepedent SNPs should behave as draws from $Binomial(2, MAF)$. Note that the standard QC procedure of filtering out genotypes that deviate from Hardy-Weinberg equilibrium is just a goodness-of-fit test against $Binomial(2, \hat{MAF})$. If you employ this filter, some LD-pruning, and still get different results, consider posting more details about the discrepancies you see - that would be indeed unexpected.

Thanks. I agree that independent SNPs behave as draws from Binomial. I think the discrepancies I observe is mainly on a population level.


text:
How to retrieve InChI key for KEGG compound?
I would like to retrieve the InChI representation for a given KEGG compound but I fail to find a direct solution for this (preferably doing it via bioservices).
One could do it via ChEBI like this:
from bioservices import *

kegg_con = KEGG()

kegg_entry = kegg_con.parse(kegg_con.get('C00047'))

chebi_con = ChEBI()

chebi_entry = chebi_con.getCompleteEntity('CHEBI:' + kegg_entry['DBLINKS']['ChEBI'])

print chebi_entry.smiles
print chebi_entry.inchi
print chebi_entry.inchiKey

which will print
NCCCC[C@H](N)C(O)=O
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1
KDXKERNSBIXSRK-YFKPBYRVSA-N

Is there a direct way which does not require to go via ChEBI?

I had a look at the KEGG rest API but couldn't find anything; I think you might need to access the data in the KEGG FTP ($$$)

@mgalardini: I checked that, but seems it is not available at all, or I overlooked it. The way above is still feasible, just requires one additional step.

It's not too bad, especially considering that you have a library to do it for you :)

I'm not aware of a functionality in the KEGG service that would allow you to do that directly. I believe your solution (using e.g. ChEBI) is the correct one.
You may have other useful service such as unichem, chemspider (you would need a login though for that one) could help as well.
Disclaimer: I'm the bioservices main author


text:
Read and write FASTA files with more information than id and sequence
I am trying to read a fasta file, manipulate is in Python (using BioPython) and then write it back.
The format of my sequences is like:
>k119_5 flag=0 multi=141.0706 len=473
AGGTTAGTCAGCACCGTTTCCGTGGTGCTGCCTTTCGCTTCAAAACCGACGGCGTCTATTACTGCATCCACGCCCCGGTGACCTGCCGTTTGTTCAATAATTGACTGTGCCGGATCGCTGTCTTCATCAAAATTAATCGGGATCGCGCCGTAGCGGTCGGCGGCGAAATGCAAGCGGTAGGGATGATGATCAACAACAAAAATCTGTTCCGCACCGAGCAACCGTGCACAGGCGATTGTCAACAATCCCACAGGACCAGCACCATAGACTGC
AACGCTTGAACCTTGTTGGATCTGCGCATTTTTTGCTGCCTGCCATGCCGTTGGCAGAATATCAGAAAGGAAAAGCGCTTTATCATCTGAAAGCAAAGGCGGTACTTTAAACGGCCCCACATTCCCTTTAGGGACGCGGACATATTCCGCCTGCCCACCAGGAACGCCGCCATACAGGTGACTATAACCAAACAATGCCGC

which is the result of the assmebly using Megahit. While trying to read, it only recognizes k119_5 as the id, while I want the whole first line as the id of the sequence. Is vital information I want to keep!
Any ideas how to work with this?

The description field in the SeqRecord object has the information you are looking for:
>> from Bio import SeqIO
>> s = SeqIO.read('genome.fasta', 'fasta') # single sequence fasta file
>> s.id
k119_5
>> s.description
k119_5 flag=0 multi=141.0706 len=473

Edit: as an aside, if you write a SeqRecord object using SeqIO's write method, the FASTA header should contain the whole original header.


text:
How to interconvert InChI and InChIKey?
I would like to retrieve IDs from several databases using InChI as an input, e.g.
InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1

One can use unichem from bioservices for this, however, these functions all require the InChIKey as input, e.g.
KDXKERNSBIXSRK-YFKPBYRVSA-N

Is it possible to interconvert the two using bioservices and if not is it possible to somehow use the functions in unichem with InChI rather than the InChIKey?
I tried:
from bioservices import *
u = UniChem()
u.get_src_compound_ids_from_inchikey('KDXKERNSBIXSRK-YFKPBYRVSA-N')

which works fine, however,
u.get_src_compound_ids_from_inchikey('InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1')

does not work and returns 400.
If it is not passible via bioservices, is there any other python API available for this?

Recently I've developed a python lib for converting InChI to InChIKey: https://github.com/liwt31/chembl_ikey. Maybe it will help.

That's a fork. Could link to the original : https://github.com/mnowotka/chembl_ikey

@gilleain it is indeed a fork. I added python3 support and fixed a bug (InChI without sminor yields wrong InChIKey).

ok, that sounds like a great contribution - but not that you developed it... Sounds like the original author could benefit from a pull request? Just a thought - I'm not the github police! :)

ChemDraw's ChemScript is a valuable tool. However, it still runs on Python 3.2. Thank you for this fork - I agree that a pull request would be very beneficial!

Not sure if directly possible in bioservices but one can do the following workaround using chemspider:
import requests

host = "http://www.chemspider.com"
getstring = "/InChI.asmx/InChIToInChIKey?inchi="
inchi = 'InChI=1S/C6H14N2O2/c7-4-2-1-3-5(8)6(9)10/h5H,1-4,7-8H2,(H,9,10)/t5-/m0/s1'

r = requests.get('{}{}{}'.format(host, getstring, inchi))
if r.ok:
    res = str(r.text.replace('<?xml version="1.0" encoding="utf-8"?>\r\n<string xmlns="http://www.chemspider.com/">', '').replace('</string>', '').strip())
else:
    print "provide a valid inchi!"

This will give the desired InChIKey
'KDXKERNSBIXSRK-YFKPBYRVSA-N'

which can be used in unichem.


text:
Filtering step for read counts data
I have around 1200 samples as columns and 60,000 genes with Htseq-Counts data. Before normalization with voom function I want to do filtering step.
I want to remove genes whose expression is == 0 in at least 10 samples.
Can I do this with read counts itself or should I transform counts to cpm?
If filtering with counts data is fine May I know how to do this.

If you have your counts in a data.frame called counts, something like this might work:
filtered.counts <- counts[rowSums(counts==0)<10, ]

For example lets assume the following data frame.
> A <- c(0,0,0,0,0)
> B <- c(0,1,0,0,0)
> C <- c(0,2,0,2,0)
> D <- c(0,5,1,1,2)
>
> counts <- data.frame(A=A, B=B, C=C, D=D)
> counts
  A B C D
1 0 0 0 0
2 0 1 2 5
3 0 0 0 1
4 0 0 2 1
5 0 0 0 2

And now you would want to remove the genes with at least 3 zeroes or more:
> filtered.counts <- counts[rowSums(counts==0)<3, ]
>
> filtered.counts
  A B C D
2 0 1 2 5
4 0 0 2 1

Every gene with 3 or more zeroes is now removed.

In at least 10 samples. So I guess it should be this filtered.counts <- counts[rowSums(counts==0)>=10, ]

No, you keep all rows with less than 10 zeroes with my example.

In your way I will be having only 1/4 th of genes remaining. I want to remove genes whose expression is == 0 in at least 10 samples. With my command now I have more than half genes remaining after filtering.

If you use it your way, you keep all the genes with more than 10 zeroes, check it for your self if you don't believe me.

So, the filtered.counts output that u gave above should be used for normalization?

If this filtering approach is what you want yes. However, in edgeR or limma manual they also give other options, such as at least 1 cpm in more than 3 samples. Usually 3 samples is then based on groups of 3 samples. If you have only 1/4 of your genes left with this filtering it might be too stringent. For limma voom, you need to have cpm values btw.

Ok. As you said I did it in edgeR way.
keep <- rowSums(cpm(counts) > 0.5) >= 10; y <- counts[keep,] So, now genes in "y" can be used for "voom" function to normalize the data. Do you think this is right?

Yes I think it is right, good luck!

10 out of 1200 samples is a very strict cut off. I'm not surprised you don't have many genes remaining.

Yes, but I'm interested in looking for a specific gene. So, to get that in my results I took that cutoff.

While this answers explains how to do it I want to address when and why and which thresholds to do it.
Filtering the genes with low counts is usually done because the counts are not reliable it would be noise, specially when there are low number of samples these genes disturb power of the analysis. However with 1200 samples having 10 or more samples without any count is valuable information because you have enough samples to estimate if this is due to a latent or known variable or if it is due to noise.
The CPM are used because the size of the library per sample affects the power of the experiment. The higher the library size the easier is to find if a gene is expressed or not. CPM are calculated to take into account that effect. Filtering by counts or CPM will result in different genes filtered if you use other thresholds than 0 (for instance I have filtered those genes that the average CPM where below 1, which would be different to filtering by those that the average counts is below 1.)


text:
True bam file encoding and viewing as text
We all have been using samtools to convert between sam and bam files.  SAM is tab-delimited file. My question resides on the ability to visualize bam files as text.  Is there any way to visualize any bam file as text in linux?
Why not work only with sam files and we need to convert them to bam?

What do you mean by "visualize bam files as text"? Do you search for something like IGV?

@Kamil S Jaron I was wondering how we can visualize as text, for example, less filename.sam would allow me to see SAM file as tab delimited.  Everyone has a good perspective on this topic.  For example @Pierre and @b.nota both mentioned that this is binary compressed file, therefore we can not see it directly without a program that interprets the bam binary encoding. An example of that program/interpreter is presented by @mgalardini and @gringer. And obviously, IGV being one of the widely used. Thank you, guys. Great answers.

Why not work only with sam files and we need to convert them to bam?

SAM is for human only, a binary format like BAM is smaller and parsed much faster by a program (and we want speed and optimize spaces for those huge files, don't we ?).
Furthermore, bam files are compressed using a  method that allow them to be accessed using random-access, that is how you can query a sub-part of bam file for a given genomic location.
you can always go back to SAM using samtools.

bam file is a compressed sam file, comparable with zip or tar. They are binary, which means we (humans) cannot read them. So you need to unpack them (convert back to sam) to read it as a text file.
bam files are smaller than sam files, which is the reason most people like to use bam files instead of sam files.

The way to visualise BAM files as text in Linux is with tools like samtools.
If you want to see the raw binary information as slightly more human-readable, you can run hexdump -C <file.bam>, but that type of text output is completely useless for further processing.
If you want more specifics about the BAM file format, you can have a look at the SAM file format specification, which also includes details on the related BAM and CRAM formats.

Assuming you want to visualize the pileup of reads across the reference, you may want to look into ASCIIgenome, a command line visualization tool for sam/bam/cram files, similar in spirit as IGV, but working in the command line.
ASCIIGenome my_sorted.bam


text:
Using the t-SNE algorithm on microarray data + an error bonus
I'm trying to use the t-SNE algorithm on some microarrays data. More specifically my data frame has 18600 columns with genes (features) and 72 rows with conditions with replicates ( 10xWt , 10xTg , etc ). The expression values are in log2 scale.
Here is the code that I'm trying to run.
# t-SNE implementation
library(Rtsne)

set.seed(1)

for(i in 1:15){

  tsne = Rtsne(data.T[,-18601], dims = 2, perplexity=i, verbose=TRUE, max_iter = 1000, pca=T)

  colors = rainbow(length(unique(data.T$classes)))
  names(colors) = unique(data.T$classes)
  plot(tsne$Y, t='n', main="tsne")
  text(tsne$Y, labels=data.T$classes, col=colors[data.T$classes])

  readline(prompt="Press [enter] to continue")
}

Please note that I'm not counting the column 18601 because this colums contains the labels/classes for each condition.
The think here is that when I execute this script, R returns me this error:

Error: protect(): protection stack overflow

Should I change the --max-pp-size or it's a bug in Rtsne package?
Also I was wondering if it is more meaningful to run the tSNE algorithm using not the log2 values of the expression level but the log fold change values in respect to the Wt (wild type) condition. I'm asking because I couldn't find a such other implementation of the tSNE on microarray data.
For the configuration of the Rtsne function I read this article
Any other suggestion on the implementation is welcomed.

Converting your data.frame to a matrix (and then removing the data.frame) will often free up enough memory that you won't run into this. Note that a matrix is more memory efficient than a data.frame and you're requiring Rtsne() to hold both in memory at the same time (many math-centric functions will end up converting things to a matrix at some point for efficiency).
For what it's worth, it's never been entirely clear to me what the interaction is between a data.frame and the pointer protection stack, but it's often the case that this solves this sort of error.

data.frames are lists of column vectors, and every one of these columns has a pointer that may need protecting. With a huge number of columns, I can easily see how this might blow the protect stack. Dense matrices are contiguous memory, they only need a single pointer.

I'd considered that, but my reading never clarified if it's really due to pointers for the columns or something else. It's unfortunate how little is mentioned about that in the documentation.


text:
How to measure euclidean distance between points with vtkDelaunay3D package?
I'm working with python vtkDelaunay3D (scipy.spatial.Delaunay) package for a special purpose. Normally the package is used for triangulation/tetrahedralisation.
How do I measure euclidean distance between points with this package?

Hi, could you provide some more details so that we can all better understand which sort of data you are working with?

I gave up using vtk package.Instead, I have decided to use scipy.spatial.delaunay which is a more robust package

@Sara If you can provide a small working example, maybe you could add your scipy-based solution as an answer.

@bli thanks but I got the answer and did that project one month ago

I think you'll find the answer here, https://stackoverflow.com/questions/20025784/how-to-visualize-3d-delaunay-triangulation-in-python


text:
How to build a BAM header file with htslib in C++?
I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly.

https://github.com/samtools/htslib/blob/develop/test/test-bcf-translate.c

The two key functions are:

bcf_hdr_append
bcf_hdr_add_sample

Now when I navigate to the header file for SAM/BAM functions:

https://github.com/samtools/htslib/blob/develop/htslib/sam.h

The only "append" function is bam_aux_append. But that's for optional fields.
My questions:

Given a bam_hdr_t header (potentially empty), how to add a new chromosome reference? (For example, add a new chrAB to the BAM file)
Given a non-empty bam_hdr_t * header, how to remove a chromosome? (For example, remove chr4 from the BAM header)

My sample code:
// Empty BAM header, nothing inside
bam_hdr_t * hdr = bam_hdr_init();

// I'd like to add a new synthetic `chrAB` into the header, how?

You're just looking in the wrong header file. The BCF stuff will be in vcf.h (you'll see the two functions you mentioned there).
With BAM files one doesn't typically need to add/remove headers, which is why there are not convenience functions for them. To add a header:

Increment hdr->ntargets
realloc(hdr->target_len, ...) and hdr->target_name.
You likely want to append the name/length to the end of hdr->target_name and hdr->target_len, respectively.

I don't think you need to update the dictionary or plain text for most uses (or at least I've never needed to).
Removing a chromosome would follow the reverse instructions, (you'd need to shift everything in the name and length arrays if you remove entries in the middle.
WARNING, doing this is prone to error!
I would strongly cautious you to be very careful when doing this. Please remember that the chromosome that entries are mapped to isn't actually stored in alignment entries themselves. Rather, they simply have an index into hdr->target_name, which saves space. Consequently, if you start adding/removing headers to/from a BAM file with alignments then you're going to have to modify all of the alignments to keep the correct alignment:chromosome associations. Otherwise you'll end up with completely useless results (this is the same caveat as with samtools reheader).


text:
How to deal with correlated SNPs in GWAS?
I have population SNP data in VCF format, and I found that some SNPs have a great similarity across samples(> 99%). For example:
CHROM_POS   s0  s1  s2  s3  s4  s5  s6  s7  s8  s9
chr1_1  A   G   G   A   G   G   A   A   G   A
chr1_2  A   G   G   A   G   G   A   A   G   T
chr1_3  C   C   C   A   C   A   C   C   C   C

The similarity between chr1_1 and chr1_2 is 0.9, because they have a different SNP in sample s9.
Is there a good way to remove these similar SNPs before I put the data into another pipeline?
Added: These SNPs will be used for GWAS analysis.

Do you mean that some SNPs have a MAF of <1%?

@Emily_Ensembl MAF is about allele frequency. But I want to remove similar SNP sequence across sampls .

The question is about if removing common SNPs a good thing.

@SmallChess, I think similar snps may will provide duplicate information and cost more resource to compute. Maybe similar snp in long distance should be kept?

It wasn't clear from your original question what you meant. Now you've edited it I can see that you're looking at haplotypic blocks with SNPs in LD. No, don't remove them. There may be causal SNPs within these blocks – you won't able to tell which SNP is causal but you can identify that the block is important.

@Emily_Ensembl, I am not familiar with haplotypic block. SNPs in this block is positional adjacent？ If they are adjacent, I would consider the distance between them, beside their similarity.

Variant that are adjacent are more likely to have alleles inherited together, due to less crossing over between them in meiosis. Pairs of variants are referred to as being in Linkage Disequilibrium (or LD) and blocks of them are called haplotypic blocks.

What you are attempting to do is known as LD-pruning.
As @Emily_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise because it would complicate interpretation: if SNP chr1_2 had some obvious function, but you removed it and are left only with its proxy chr1_1 (which could be far away), you will have a hard time connecting the dots. Typical GWAS tests are handled pretty quickly on modern computers, and are linear in SNPs, so there is no computational reason to prune either.
On the other hand, many fancier procedures would require LD-pruned data, because they either
a) are very computationally expensive, or
b) assume independent SNPs.
I suggest reading PLINK's LD pruning documentation, as it pretty much covers the standard ways to do this.


text:
genes with highest RPKM/FPKM for a human RNA-seq experiment?
Which genes/transcripts show up as having the highest RPKM/FPKM on a human RNA-seq experiment? Is it usually the same genes/transcripts?
EDIT: mRNAs (e.g. poly-A RNA-seq preps) and not including small RNAs.

What effort did you take? Did you look at available data from GEO already?

I'm guessing there's some missing context here of the form, "We've done a bunch of RNAseq experiments and X, Y, Z always show up as having the highest FPKMs. Have other people seen this or are we just doing something wrong?" If so, please include it.

@DevonRyan the context is more: "here is a handful of genes that always come up on top for a human sample. If you don't see them, something may have happened".

In that case, you may be looking for housekeeping genes -- this doesn't quite answer the question as stated (they aren't always the highest expressed, but they are usually present).

As asked as it is, the answer is probably no.
Indeed, the highest expressed RPKM/FPKMs will be different from one condition (or tissue) to another one for example.
Then, you may also have technical artefacts due to the wet-lab part or to the normalisation. For example, mitochondrial genes are often reported in the top expressed genes.
Now, if you want to compare your results to an outside source, I recommend checking the EBI gene expression Atlas* and more specifically the Baseline expression experiments.
Pro-tips:

You can filter the genes based on their nature (eg, mRNA) and/or their expression.
You can also visualise the overall expression of a specific gene across several experiments.

*Disclaimer: I am related to this project myself.


text:
Missing data mappings in mygene.info while trying to convert Genes Ensembl Ids to Entrez Ids
I need to convert a lot of Ensembl Ids to the relative counterpart in Entrez (e.g., ENSG00000157764 > 673).
I found mygene.info and it seems what I needed. Let's see the query about ENSG00000157764 in action. We can easily find the key-value "entrezgene": 673
Unfortunately, the mapping seems not complete in some cases or, maybe, I'm missing something...
Example:

Gene: NDUFV2P1 (NADH:ubiquinone oxidoreductase core subunit V2
pseudogene 1)
Expected Mapping
ENSG00000267809 > 4730

If I query mygene for ENSG00000267809, I get this result, but this time no data about the entrez id.
What it sounds pretty odd is that if I query another mygene's endpoint for NDUFV2P1, 2 identical records are returned, respectively with Ensembl and Entrez Ids.
Is there a way to leverage such a service to get clean data? Otherwise, is there another robust conversion mapping strategy out there?

"I found mygene.info and it seems what I needed" : you should try biomart http://www.ensembl.org/biomart/martview/ .

Cannot mygene be queried for my use case? Biomart seems very user unfriendly from the API point of view

mygene is a no-name site of uncertain quality. Biomart is a definitive source for reference information. A mild incovenvience in API shouldn't be a deciding factor (especially when you can just download a small text file with everything).

ok @DevonRyan thanks for your suggestion. I'll try to find out the text file with the mapping that I need.

As others have mentioned in their answer, bioMart is usually the best place to go for this infomation because it draws its data directly from the Ensembl database.
However, you will find you do not get a full translation of Ensembl IDs to Entrez IDs even there. The reason for this is there simply does not exist a simple 1-to-1 mapping of Entrez IDs to Ensembl IDs: they are not just different names for the same thing, but are actually different things. At some point someone (or more likely some algo) has decided that the structure Ensembl ID A is referring to is similar enough to the structure that Enztrez ID B is referring to to link the two together, but such a linking is not always possible.
In the example you gave, the Ensembl database does not list ENSG00000267809 as being equivalent to 4730, neither does the Entrez database list the opposite, even though both genes have the same name. Because of this, it is very unlikely that any service you use to convert the IDs will pick up this match.
In this particular case, the reason for this is probably that the gene is a processed pseduo-gene and is not annotated in Refseq/NCBI as producing any transcripts.

Thank you once again Ian. Very clear. Therefore, If I understand well, the problem I've reported in the question is not real. I mean, mygene is doing a good job keeping those two concepts as separate items.


text:
Adjusting phenotypes by regressing out covariates
I'm trying to use the bfGWAS tool, which analyses GWAS data and integrates functional annotations to identify casual SNPs (paper and
github). In the user manual, it states:

We recommend first regressing covariants out from the original
  phenotypes, and then provide bfGWAS the corrected phenotypes (i.e.,
  residuals from the regression model with covariates).

It is unclear to me what this means and how to go about determining the corrected phenotypes. Typically, in an association test with covariates, a regression would be performed using a model like
pheno ~ snp + covariate

I could then get the residuals from this model. However, I'm not certain how to go about it for all SNPs, and then adjust the phenotypes. I have case-control data, so the phenotypes are binary.

Instead of the "observed phenotype", it's asking for a "corrected phenotype", which is the residuals from pheno ~ covariates. In R, one would get that as follows:
fit = lm(phenotype ~ covariates, data)
corrected_phenotype = resid(fit)

Of course covariates is just a place hold for something like age + gender or something along those lines.

Hi @Devon Ryan, I have a question regarding phenotype. Will you be able to help. only if you have time. Else you can ignore it. https://bioinformatics.stackexchange.com/questions/11028/use-of-electronic-phenotype-in-ehr


text:
How to determine the primary Uniprot accession number from a set of accession numbers?
Given a list of Uniprot IDs that are linked to an Ensembl gene ID, is there a way to systematically determine the UniProt accession corresponding to the protein product of the gene's primary transcript with no other information?
According to ExPasy

Researchers who wish to cite entries in their publications should
  always cite the first accession number. This is commonly referred to
  as the 'primary accession number'. 'Secondary accession numbers' are
  sorted alphanumerically.

But what if the order has been jumbled or compiled from a different source and resorted?
For example:
Ensembl:
ENSMUSG00000035642
Uniprot:
Q8R0P4, Q8CF11, D6RJK8, D6RJJ4, D3Z442, D3Z1Q3, D3YZD8, D3YY39, D3YX09, D3YWY5
This question is cross listed on the Biology stack-exchange.

If I understand this correctly, each Uniprot entry has the same primary accession number (which you can find in the "Entry Information" tab of the Uniprot entry; ctrl+f "primary"). Some of the IDs that you give here describe different sequences (which thus have different Uniprot entries and different primary accession numbers). I have checked the first 5 IDs you show here. They all map to different sequences, and should have different accession numbers. They are isoforms, mutated or truncated sequences and God knows what else, but they are definitely not the same protein sequence.

I never claimed these are representative of the same protein squence. The question above states: Given a list of Uniprot IDs that are linked to an Ensembl ID

No, you didn't :) What I am trying to say is that the way I understand it there shouldn't be an 1:1 mapping between Ensembl and UniProt in the first place. They describe different things. Am I missing something? It seems to me like you are trying to pick one protein ID to connect to one gene ID.

@galicae didn't do what? I'm trying to determine if there is a simple procedure to identify the primary Uniprot accession number from a list of associated accession numbers with no other outside information. I don't think I can state it more clearly.

@JamesDraper do you mean the primary UniProt accession of the protein product of the primary transcript associated with that gene? Or do you mean the primary UniProt accession of each of the protein isoforms encoded by the gene? Please [edit] your question and clarify what you need since there is often (usually) a one-to-many relationship between genes and proteins.

@JamesDraper you said you "never claimed these are representative". I agreed :)

@terdon the first option you mentioned. Feel free to edit the post as you see fit.

@galicae I'm glad that we are on the same page.

@JamesDraper OK, and how do you need to do this? Manually? Does it need to be scalable to hundreds of accessions? You mention "no other information", does that mean that going from the Ensembl Gene ID, to the Ensembl Transcript ID for the primary transcript and then getting the UniProt accession(s) associated with that is not an option?

@terdon honestly I was was hoping that this question had a simple obvious answer that I just didn't know like: sort in in descending order and take the first member of that list. Or if the something like the primary accession could be recognized by  some kind of special regex. However it has become abundantly clear that the answer to the question is no there is no way to determine the primary accession from the a list of associated accession numbers with no outside info. If anyone would like to restate that as answer I'll accept and upvote.

Reposting my answer from the related ticket in the Biology section:
I think there is an issue with the terminology. The "primary" accession number, is the first accession number in cases where an entry has more than one accession number, as described in http://www.uniprot.org/help/accession_numbers:

Entries can have more than one accession number. This can be due to
  two distinct mechanisms:
a) When two or more entries are merged, the accession numbers from all entries are kept. The first accession number is referred to as the

‘Primary (citable) accession number’, while the others are referred to
  as ‘Secondary accession numbers’. These are listed in alphanumerical
  order.
b) If an existing entry is split into two or more entries (‘demerged’), new ‘primary’ accession numbers are attributed to all

the split entries while all original accession numbers are retained as
  ‘secondary’ accession numbers.
Example: P29358 which has been ‘demerged’ into P68250 and P68251.

Both reviewed and unreviewed entries can have primary accession numbers.
What you probably mean, as previous posters understood, are accession numbers of reviewed entries as opposed to unreviewed ones.
In that case, you can indeed add "reviewed:yes" to your query, e.g. when you are using the UniProt ID mapping, http://www.uniprot.org/help/uploadlists

Consider the circumstances in which entries will have more than one accession (taken from the same page you linked in the original post).

Entries will have more than one accession number if they have been merged or split. For example, when two entries are merged into one, the accession numbers from both entries are stored in the AC line(s).
If an existing entry is split into two or more entries (a rare occurrence), the original accession numbers are retained in all the derived entries and a new primary accession number is added to all the entries.

As mentioned in the comments, this does not imply that there is or should be a single "primary" UniProt accession associated with (for example) an Ensembl gene ID. In some cases, a single protein record may have multiple accessions due to annotation changes, and the suggestion from your original post is to only report the primary accession. However, there may be multiple proteins associated with an Ensembl gene ID, and each protein will have its own accession (or set of accessions).
Choosing a single representative isoform for a gene is an entirely different question.

Not with no other information, no. There is no way of looking at a list of accessions and guessing which one is primary. Primary accessions are just accessions like any other; they don't have any specific characteristics.
This is further complicated by the fact that each gene can (and usually does, for higher eukaryotes) have multiple transcripts, each of which can produce a different protein product. So in order to go from gene identifier to a protein, you need to first find the relevant transcript.
So, as a first step, get the primary transcript IDs for each of your genes and then retrieve the protein accessions associated to them.


text:
Interpreting Intergrative Genomic Viewer (IGV)
I was following a tutorial on "Tuxedo Genome Guided Transcriptome Assembly Workshop" and was wondering how to interpret the following:

From what I understand from 'Color Legends', the color blue represents something that is below normal and red represents something that is above normal.
Also, my understanding is that merged.gtf is the reconstructed transcript and genes.bed is the reference. They look to be similar for the most part, but seems to have breaks between. What's the appropriate way of interpreting this result?

You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad institute software engineers replying here.
Your next question is likely why the color of the reads and the orientation of the genes are opposite of each other. This is because in modern directional protocols, read #2 in a pair matches the orientation of the gene from which it arose.

Thank you for your reply. I'm not sure I understand what you're referring to by read #1 and read #2 based on the picture I attached.

Ignore that part if you don't understand it. The direct answer to your question is that the colors represent the strand from which the reads arose.


text:
Where to download JASPAR TFBS motif bed file?
I am interested in determining if any transcription factor binding site motifs are enriched in some BED files from a DNA methylation experiment.
I am looking for a database that has BED Files containing regions enriched for specific transcription factor binding motifs in homo sapiens.
I know that JASPAR contains TFBS motifs, but I can't find a download for regions of hg19 that are enriched with those motifs.
EDIT:  I suppose I could determine TF motif enriched regions using HOMER if there are no curated databases with this information, but I would rather not.  Your help is appreciated!

You could make your own such file with FIMO, the JASPAR MEME files available via the MEME download site, and your genome build of interest (per-chromosome FASTA files from UCSC goldenpath, say), e.g., for assembly hg19 and UCSC chromosome naming scheme:
$ for chr in `seq 1 22` X Y; do echo ${chr}; wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg19/chromosomes/chr${chr}.fa.gz | gunzip -c - > hg19.chr${chr}.fa; done

For each chromosome (omitting the for loop code for brevity), you can use the fasta-get-markov in the MEME suite to generate a 5th-order Markov model of kmer frequencies (why 5th order?):
$ fasta-get-markov -m 5 hg19.chr${chr}.fa > hg19.chr${chr}.5th_order_background_model.txt

Here is one way to download a non-redundant, vertebrate JASPAR database (MEME-formatted):
$ wget http://jaspar.genereg.net/download/CORE/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.txt

To run a FIMO query with the 5th-order background model and JASPAR database, using a p-value threshold of 1e-4:
$ fimo --verbosity 1 --bgfile hg19.chr${chr}.5th_order_background_model.txt --thresh 1e-4 --text JASPAR2018_CORE_vertebrates_non-redundant_pfms_meme.txt hg19.chr${chr}.fa \
  | tail -n+2 \
  | awk -vOFS="\t" -vFS="\t" '{ print $3, $4, $5, $1":"$2">"$10, $8, $6 }' \
  | sort-bed - \
  | starch --omit-signature - \
  > hg19.chr${chr}.bed.starch

The output of a FIMO search is a text file which can be modified into a sorted BED file via sort-bed, an example of which is shown above.
To get a faster answer, each per-chromosome query can be put onto its own node on a computational cluster.
Using a p-value threshold of 1e-4 will result in very large files (when uncompressed). Using tools to compress BED files is useful here, as is using a more stringent threshold (1e-5, etc.).
Use of BEDOPS starch to compress BED files (as shown above, for instance) will give better compression performance than bgzip with per-chromosome random access and integration with BEDOPS bedops and bedmap CLI tools. The bgzip option offers integration with tabix and offers finer-grained random access, but it will use more disk space and doesn't integrate with BEDOPS tools.

This is a decent answer, so I am accepting it.  If you ever hear of these files being added to a database or of a database of this type of data please let me know :)

I can't presume to speak for groups like UCSC or Ensembl, but I'd guess that each analysis of this kind is short-lived, based on the fact that it really comes down to what motif models are used (which get updated), what statistical parameters get used (which are down to the experimenter), and the reference genomes involved. It's perhaps not something unique enough that you'd generate once in a blue moon, that it could be efficiently made available in a genome browser or via a biomart-like interface. It is basically something you make to meet your particular experimental needs, I think.


text:
Converting Ensembl Gene IDs to Entrez Gene IDs through biomart
Well, I'm trying to convert a list of Human Gene referenced by Ensembl Gene IDs to Entrez Gene IDs. I have been advised to use biomart.
I tried to get a kind of conversion table for all human genes. I don't know if my settings are wrong, but I didn't find any checkbox for Entrez Gene Id in the section Attributes > External References. I've just found the Entrez transcript name ID (see below), but it's not what I need.

How can I use biomart for such a conversion?
Edit:
According to this video tutorial EntrezGene ID should be an option in biomart (see screenshot below).
That video tutorial is outdated.

That tutorial is 8 years old.

Yes, @Emily_Ensembl. I'm sorry, I had gotten it from an old documentation page : )

The ID you need is the NCBI gene ID, which is the same as the EntrezGene ID.

Thank you Ian, it seems working. Anyway, there is no a one-to-one matching (i.e., some Ensembl IDs do not have an NCBI counterpart). Is it normal?

Yes. See my answer to your other question. https://bioinformatics.stackexchange.com/a/2474/235

Conversion using R:
library(biomaRt)

mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))

    genes <- getBM(
      filters="ensembl_gene_id",
      attributes=c("ensembl_gene_id", "entrezgene_id"),
      values=ensembl.genes,
      mart=mart)

Where ensembl.genes is a vector of Ensembl gene IDs.

As of 2021-08-05, the attributes argument should be set to attributes = c("ensembl_gene_id", "entrezgene_id"). I believe that the "_id" appendage has been added to the attribute since this answer was first posted.


text:
How to represent a deletion at position 1 in a VCF file?
I am writing a small script to write a VCF file from short multiple sequence alignments, and stumbled into an interesting corner case.
This is the way a deletion should be represented in VCF format (I assume version 4+):
NC_000016.9  2138199 rs1799758   GTGAG   G  .   .  RS=1799758;...

But what happens if the deletion is at the first position? Immagine the following two aligned sequences:
>reference
AGTC
>sample
-GTC

How should it be encoded in VCF format? According to this thread on SEQAnswers, it might be necessary to include the first position after the end of the deletion, like this:
chromosome  1  AG G  .  .  [...]

Is this a VCF 4+ correct way to encode such variants?

From the latest VCF specs (page 8):

REF - reference base(s): Each base must be one of A,C,G,T,N (case insensitive). Multiple bases are permitted. The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event occurs at position 1 on the contig in which case it must include the base after the event;

So I think you're right


text:
Extract LINE-1 methylation status from WGBS data
I've been trying to figure out the correct way of extracting methylation status of particular elements from WGBS data, and found several possible approaches, but would like to discuss which one could be the best:

using software dedicated to the extraction of these elements on WGS (like TraFiC or MELT and then quantify methylation (like Bismark)
using bedtools to scan bam file with a filtered version of UCSC RepeatMasker table in bed format) and then quantify methylation
I've seen the use of custom perl scripts on fastq files like in this paper, in which the Supplemental Methods Section states: "A custom Perl script was written to assay methylation of CpGs within a previously published motif of LINE-1 elements (Yang et al., 2004). Methylation levels for each of the four CpGs within the consensus LINE-1 element were determined globally for each WGBS sample prior to alignment using Line1_FASTQ.pl. Dup15q and control were compared and tested for significance using an unpaired t-test."

Any thoughts on this?

Have you checked out WGBS Tools on github? All of the custom scripts used in the paper that you referenced are in that repository. Each tool has a brief description of how to use it.

After a more extensive search, I found a useful answer that I'll share in case anyone else wonders how to do this. I'm using MethylKit package for differential methylation expression and it is quite easy to integrate with Genomic Ranges and rtracklayer in R to extract particular element methylation differences:
element.of.interest <- import.bed("FILENAME.bed") #download custom track from Genome Browser
element.meth <- diffmeth[as(diffmeth, "GRanges") %over% element.of.interest]

where "diffmeth" is the resulting MethylKit object from the differential methylation analysis.

Via BEDOPS:
$ bedmap —-echo —-echo-map-id-uniq <(bam2bed < reads.bam) <(rmsk2bed < repeats.rmsk) > answer.bed


text:
Transform traditional blast output to `--outfmt 6`
I have run a blastx of metagenomic databases (raw illumina reads) using the nr database. Unfortunately, I forgot to add the --outfmt 6 argument to the code and got the traditional output.
Could I parse this output into --outfmt 6?

By "metagenomic databases" do you mean collections of reads or assembled contigs? BLAST searches of metagenomic reads vs nr is a very bad idea. Much better tools exist for classifying, etc.

@ Daniel Standage  I refer to a collection of metagenomic reads. Could you suggest any other classification method for shotgun metagenomic reads?

Sure! Although you should take my recommendations with a grain of salt, since this isn't really my area of expertise. But tools like Kraken and Kaiju seem to be pretty popular.

I think the standard BLAST output is a moving target so the exact parsing will depend on the version of BLAST you use; anyway I would recommend re-running with the --outfmt 6 flag unless you don't have the computational resources. If it's too slow you can use DIAMOND as a replacement for BLAST with very similar commands

Biopython has a parser for Blast's default text format (which I assume is the one you have). Just use the blast-text string to specify this format in SearchIO's parse method.
from Bio import SearchIO
results = SearchIO.parse('myfile.txt', 'blast-text')
for entry in results:
    # extract the information you want
    query = entry.id
    for hit in entry.hits:
        subject = hit.id
        for hsp in hit.hsps:
            identity = hsp.ident_num *100 / entry.seq_len
            aln_len = hsp.aln.get_alignment_length()
            mismatches = aln_len - hsp.ident_num # not sure about this
            gaps = hsp.gap_num
            q_start = hsp.query_start + 1
            q_end = hsp.query_end
            s_start = hsp.hit_start + 1
            s_end = hsp.hit_end
            e_value = hsp.evalue
            bit_score = int(hsp.bitscore)
            # print in outfmt 6 format
            print('\t'.join([str(x) for x in [query,
                                      subject,
                                      identity,
                                      aln_len,
                                      mismatches,
                                      gaps,
                                      q_start,
                                      q_end,
                                      s_start,
                                      s_end,
                                      e_value,
                                      bit_score]]))

Have a look at Biopython's docs for more information on how to parse Blast results.

This is a very versatile approach, thanks.

@andresito, please remember to accept the answer that helped you the most

The blast_formatter command allows you to produce output in any of the formats supported by BLAST. However, it requires that you run blastx with --outfmt 11 (ASN.1 format). I'm not sure there are any tools that will do the conversion from the full default BLAST report. This could potentially be scripted, but would be complex enough that you're probably better off re-running BLAST.


text:
Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?
For the purpose of testing phylogenetics software, we're using synthetic sequence data (i.e., it's computer generated, and there's no actual biological organism with that e.g. DNA sequence).  For real-world biological sequence data, we would say the sequence belongs to a "species" or "taxon" (or "taxonomic group"), but this feels inaccurate since there is no actual species or taxon.
Question: Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?
I'm wondering if this is standard practice, despite the synthetic nature of the data.

How were the synthetic reads generated exactly?

Using evolver, packaged with PAML.

Phylogenies can be made with/of non-species groups. I think "taxon" is perfectly fine. You might also consider the use of the "leaf / leaves" or "terminal" if you want to emphasise the tip nature of the organisms. Even "sequence" would be fine.
Ninja edit to my own suggestion: since you're making a phylogeny from sequences of unknown context or relationship (i.e. you can't say that they're sequences from a "species" or "OTU" or from anywhere), agnostic terms like "taxon" or "sequence" are more correct.

'Taxon' seems to be frequently used in the literature for simulated sequences. See this or this for example

When looking at genomic data for bacteria, taxa above the species level are operational taxonomic units (OTUs). OTUs are designated based on percent identity of the sequence (specifically the 16S rRNA gene sequences in most cases). The classification is usually based on rule of thumb, and some rough guidelines are:

97% for Genus
94% for Family
88% for Order

Using OTUs is common practice. I don't see much of an ideological difference using OTUs with synthetic data.
Konstantinidis, K. T., & Tiedje, J. M. (2005). Towards a Genome-Based Taxonomy for Prokaryotes. Journal of Bacteriology, 187(18), 6258–6264. http://doi.org/10.1128/JB.187.18.6258-6264.2005

OTUs are restricted to above the species level: they can be used for any group of closely related individuals. The use of OTU is meant to be agnostic as to actual group or relationship, underlining that the grouping was arrived at "operationally", i.e. methodologically. In any event, my feeling is that the context in which OTU is usually used is fairly distinct to a generic phylogenetic problem.


text:
VCF merge containing CNV
How do I merge VCFs files containing CNVs?
I use vcf-merge, a VCFtools function, and after bgzip and tabix, SAMtools, to index and tab separate variants, but I don't know if it is the right way.

Can you post an example line from the VCF? If the CNV is formatted like a SNP, just with insertion/deletion alleles, I see no reason why there should be any problem. Then again, you can always just extract a line from the merged file and see if it looks as expected.

vcf-merge (or better: bcftools merge) does a column-based merge of VCF files, basically taking on an additional genotype column to each variant. While this sort-of works, it doesn't take into account all the information about variation. There might be low coverage at a region in one VCF file that indicates a deletion, but the deletion is only observable when considering other samples with high coverage at that same region. It's also possible that INDELs could be defined differently in different samples, requiring an additional bcftools norm step for left-normalisation prior to the merge.
Assuming this column-based merge is what you want (and if not, then bcftools concat is your friend for row-based merging), then it's better to regenerate the merged VCF file as a multi-column file directly. One way of doing this is feeding all BAM files into samtools mpileup to generate genotype likelihoods, then feeding the result of that into bcftools call. Something like the following:
samtools mpileup -uf ref.fa in1.bam in2.bam in3.bam | \
  bcftools call -mv -O z -o variants_ref_1-3.vcf.gz

depending on your goals you could possibly use tools like https://github.com/mkirsche/Jasmine or simple bedtools intersect on the vcf entries could help (e.g. harmonize after the calling)
or if you have a cohort maybe possibly doing joint-CNV calling or "genotyping CNV" in a cohort is useful http://penncnv.openbioinformatics.org/en/latest/user-guide/joint/ or https://github.com/abyzovlab/CNVpytor#merging-calls-form-multiple-samples


text:
Server for finding kmers in set of sequences
Is there a server/website somewhere where I can submit a list of DNA/RNA sequences and find the list of kmers hits and organisms where it's found? I checked the kraken website, but they don't have a webserver for it.
I am specifically looking for kmers, as a complement of just doing BLAST on the NCBI nt database or similar.
EDIT: So far I have seen:

the Kraken app on Illumina BaseSpace, but it's not trivial to upload FASTQ data that hasn't been produced from an Illumina instrument.
the CLARK trio of apps in InsideDNA, tried once with default parameters but produced empty results.

Why you need a webservice? Would not be easier to calculate them by your own?

One Codex is powerful and user friendly, and works in a similar way to Kraken underneath
https://www.onecodex.com/
The Kaiju web server might be of interest and is more sensitive than exclusively k-mer based methods
http://kaiju.binf.ku.dk/server
The title of this question is misleading – it sounds like you're actually looking for a webserver for taxonomic classification using k-mers.

Ah yes! I now remember having used onecodex in the past. I'll try both.

I tried onecodex both uploading a file and also copy+pasting the sequence into the 'Search' menu, and both came with no matches or 'Error'. Now trying kaiju...

I tried kaiju with standard parameters, no hits. I tried with relaxed parameters and got some hits. It did what it needed to do, and quite quickly as well. Thanks!

Glad that worked for you. k-mer based approaches have limited sensitivity in the presence of genetic variation and/or sequencing errors. Kaiju isn't exclusively k-mer based unlike One Codex IIRC. You could always try a more sensitive approach such as BLAST.

Perhaps you could edit the title to make this question more useful to others.


text:
Using BLAT command line tool to blat split sequences
I have the nucleotide sequence:
AATTGAGGCACATTTTTTTTTAGACAGTCTTGCTCTGTTGCCCAGGCTGGAGTGCAGTGGTGTGATCATAGCTCACTGCAGCCTCGACCTCCTGGGCTCAACAAAGCACACAGTGGGCGGATCCCCACCAG
When I blat this on UCSC Genome Browser, the first hit is a match which spans 6572 basepairs on chromosome 19, and essentially matches the first portion of the sequence to one spot on the genome, and the second portion to another part of the genome approximately 6.5kb downstream.
However when I use the BLAT command line tool, it doesn't find the second half of this match, it only matches the first portion.
Here are the flags I'm using:
-t=dna
-q=dna
-minIdentity=90
-out=blast8
How do I get the command line tool to mimic the behavior of BLAT on UCSC Genome Browser?

The FAQ offers an answer:

I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.
We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. The web-based blat is tuned permissively with a minimum cut-off score of 20, which will display most of the alignments. Other than to confirm that your command-line blat is working, there is little use in perfectly replicating the web-based blat results. We advise deciding which filtering parameters make the most sense for the experiment or analysis. Often these settings will be different and more stringent than those of the web-based blat. With that in mind, use the following settings to replicate the search results of the web-based blat:
[ . . . ]
standalone blat:
blat search:
blat -stepSize=5 -repMatch=2253 -minScore=0 -minIdentity=0 database.2bit query.fa output.psl


text:
Clarification on Gene Enrichment
When I run a GSEA analysis on two conditions from the same RNaseq (negative control PBS injection VS positive control CpG injection) from the same dataset/same gene list, I get results that look something like this:

Notice in my example that many gene sets are significantly enriched in PBS VS CpG, but none are very significantly enriched in the inverse, CpG VS PBS.
I have a basic question:
If we're comparing two items, say A and B, shouldn't a up-regulated gene set that's enriched in A have a corresponding down-regulated gene set that's enriched in B?
I'm consistently confused about how to interpret the fact that one of my two conditions has many more enriched gene sets than the other. What am I missing/misunderstanding?
Thanks!

If you compare A vs B the genes's fold change will have the opposite sign to B vs A. So will be the gene set up or down-regulated depending on the comparison to take.
The gene set test analyze if a given group of elements is sorted in a certain way in a list (I am talking about the GSE like the one you performed or the one in the Broad Institute). Usually it is used to say that gene set X is up-regulated in the comparison A vs B, that is, gene set X (in general) is more expressed in A than in B. The enrichment is measured in a enrichment score (or a normalized enrichment score, the higher, the clearer trend/distribution group X has. To asses X is really differentially expressed a p-value is also calculated.
Despite that the number of the up-regulated gene sets sum up all the gene sets tested as pointed by juod, the number of gene sets with a a p-value below 5% and 1% is not the same in both comparisons which indicates that you are not really using the same samples in both gene sets.
Check how you performed your comparison. Are the genes in each comparison with the same value and opposite sign?
Check also the software used to run the GSEA analysis (if I recongize the type of report is from a package in Bioconductor [I don't recognize which one now]) as you (or me) might have misunderstood something or there is an error if the comparisons are done properly. For instance it could be that the report of p-values below the threshold would be from only the up-regulated gene sets.

Thanks for the answer! I'm still having some trouble, most likely because I'm so new to this. This came from RNAseq data from the same experiment, but there are 3 "PBS-Treated" samples and 3 "CpG-Treated" samples. I'd expect lots of upregulation in the CpG treated samples.
In the GSE results, one of the top PBS groups is 'GO_NEGATIVE_REGULATION_OF_B_CELL_ACTIVATION'. This means that particular gene set is enriched in PBS and not in CpG, right? Wouldn't I expect to see B cell activation in the CpG sample, then?

It might make it easier to explain the biological answer I'm looking for. I'd want to be able to say, "When we injected with CpG (as opposed to PBS control), the following gene sets were upregulated:". Is it just that there are no such gene sets?

@JulianStanley If you use GO as gene sets you don't use the relationship and the structure of the gene ontologies, (use topGO Bioconductor packages that take that into account). I think the GSE is not correctly performed: did you check what I said? But yes GSE will answer this kind of questions when done correctly. See my edit too

@Llopis I think the key here is to find out which test exactly is performed in that package. Your description sounds like Mann-Whitney or something similar, which I agree should be symmetric. However, I am more used to enrichment testing by Fisher's test (i.e. contingency table "belongs to set-does not" vs. "is upregulated-is not"). Maybe the latter approach doesn't necessarily give symmetric p-values?..

@juod I agree that knowing which test OP used is important. The Mann-Whitney is different from GSEA in that you don't have up/down-regulated gene sets, but over represented (OR) gene sets. The p-values if it's a OR test won't be symmetric, if it is a GSEA they should.

Hey all--thanks again for the discussion. I just use GSEA (the software provided by the Broad). With basically default settings--including noise-to-signal stats and I used the GO gene sets from MsigDB. You also said to check that the values are the same but opposite, but I wasn't sure what you meant by that. They're RNAseq counts--so no negative values.

@JulianStanley Great to know the tool. Did you provide your expression or a sorted list of differentially expressed genes? I said to check the opposite sign thinking you have a fold change or as it is called in the program a sorted list. If you upload the AvsB and the BvsA the same gene set will have the same p-value and different sign, 0.5 to -0.5 normalized enrichment score for instance.

@Llopis ahh that makes sense! The program did provide me a gene list with enrichment scores, but I didn't know how to interpret it. I just provided my expression (deseq normalized). Notably I also tried AvB and BvA and the results looked the same at first glance

@JulianStanley Great to know that finally it make sense, I'll try to update the answer.  If you have other questions I'll be glad to answer them. Also to mark as solved consider marking one answer as accepted.

While I'm not familiar with GSEA software in particular, I believe your problem is that it only tests for upregulated gene sets. Notice:

866/4408 gene sets are upregulated
3542/4408 gene sets are upregulated

3542+866 = 4408. I.e., 866 sets have higher mean expression in positive condition, the rest have higher mean expression in negative condition. To avoid this confusion, I personally would reserve using the terms "up/down-regulated" only for sets where this change is at least remotely significant.


text:
SNP-phenotype association analysis only using the SNPs on a specific gene
I know that the GWAS association p-value threshold is 1e-8. This makes sense because in GWAS, you make ~1 million hypothesis (i.e. use that many SNPs in the association test). However, let's say I do not do genome-wide analysis but only use the 100 SNPs from a specific gene to check for association with my phenotype. Then, do you think a raw p-value of 1e-4 (which becomes 1e-2 after multiple hypothesis correction) is a reliable p-value to come to a conclusion that there is an association between that SNP and the checked phenotype? Obviously, if I checked all 1 million SNPs instead, this p-value would be insignificant after correction, but I am not checking 1 million SNPs, as I mentioned.

I know that the GWAS association p-value threshold is 1e-8

This may be a common threshold of statistical significance that is used, but it's definitely not an absolute value. It's a hack to try to work around many issues with GWAS associated with testing millions of SNPs. Unfortunately, the most relevant issue in GWAS (for spurious significance) is unexpected shared genetic structure in the cases or controls, and this cannot be excluded by a p-value threshold, regardless of how low it is set.
In your 100-SNP example, you are correct that the bonferroni correction would be a threshold of 1e-4 (for a desired significance threshold of 0.01), but it's important to be cautious about GWAS results, even when the p-value indicates otherwise. Trust no one, especially yourself. P-values should never be used as an indication of importance (which suggests that most GWAS results should probably be reconsidered), and are better used in conjunction with other evidence for identifying relevant / significant tests.
My recommendation for GWAS (or any similar bulk genotyping study, as in your 100-SNP subset) is that at least a bootstrap sub-sampling process is carried out, to make sure that observations found within the comparison of interest are at least well replicated when comparing smaller sub-groups of the cases and controls. More information on that can be found in the poster I presented at Queenstown Research Week last year, and detail about the implementation in my draft preprint (which will probably never be properly published, given how old the research is).

Thanks for the reply, those make a lot of sense. A followup question: I was also told by a GWAS expert that when the minor allele frequency is low, one would not trust the p-value even if it is significant. He said so when I showed him the p-value of 1e-4 for a SNP with a 2% minor allele frequency. Do you think it is due to the need for the additional bootstrap analysis you suggested? I guess a low p-value would not be replicated in bootstrap analysis in case of rare variants and small overall sample size. Could you share your suggestions specifically for analyzing rare variants?

Whether or not it will be internally replicable will depend on the sample size of the subsampled population. It would be possible to specifically subsample to increase the frequency of single SNPs, but that might cause problems with associated interaction effects.

Certainly, you need to take more care about rare variants. Most tests will provide both a p-value and an effect size (which is better than the p-value for ranking SNPs). If the effect size is large, then it lends more weight to an association near the SNP locus. It may also be a good idea to do a haplotype association test within the haplotype block(s) that include the SNP, as any linked SNP could be causative. Bear in mind that GWAS are only testing for association, and there are other follow-up studies that should be done to support findings of association.


text:
How to get a list of genes corresponding to the list of SNPs (rs ids)?
Is there a way for me to get a list of genes given a list of SNP rs ids? I found several questions asked with a similar goal years ago, and the answers are always about using multiple online tools, and/or R programming.
I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )
The relevant question I found (that are asked and answered more than five years ago.)

https://support.bioconductor.org/p/39278/
https://www.biostars.org/p/10437/
https://www.biostars.org/p/413/

Actually functional annotation is a key step in SNP calling pipelines. I have used SnpEff and it worked perfectly. Even if you have to build your own data base given the reference genome and the gene models against the SNPs were called. I have a script if you want to use it as guide.

I wonder if there is a simpler solution recently? (and hopefully, I can solve it within the scope of python. )

A simpler solution, I don't know... but this is at least one Python solution using Biopython's ELink method via NCBI's Entrez E-utils.
The Biopython library is flexible enough (they have an in-depth tutorial worth reading) to modify the code below to to fit project-specific needs.  For example, the mapping between snp id and gene id is lost with this bulk request method.
from Bio import Entrez
Entrez.email = "youremail@gmail.com"

snp_ids = ["rs8179065",
           "rs12899845",
           "rs8179066",
           "rs80357332",
           "rs80356920"]

record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=",".join(snp_ids).replace('rs', ''),
                                  db="gene"))

This returns a record object
[{'LinkSetDb': [{'DbTo': 'gene', 'LinkName': 'snp_gene', 'Link': [{'Id': '374654'}, {'Id': '90381'}, {'Id': '672'}]}], 'ERROR': [], 'LinkSetDbHistory': [], 'DbFrom': 'snp', 'IdList': ['8179065', '8179066', '12899845', '80356920', '80357332']}]

which contains several iterable fields -- 'Link' has a list of the Gene UID's that can be converted to names using ESummary.
for gene_id in record[0]['LinkSetDb'][0]['Link']:
    handle = Entrez.esummary(db="gene", id=gene_id['Id'])
    uid_record = Entrez.read(handle)
    handle.close()

    uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
    print (gene_id['Id'], uid_summary['Name'])
    # ...or something more interesting than print

374654 KIF7
90381 TICRR
672 BRCA1

An option to preserve discrete mapping, is to loop over each snp id and extract the gene(s):
although I don't endorse the method below, especially for long lists of snps, you may flood the server with too many requests, maybe someone more familiar can comment on request limits or best practices for Entrez
for snp_id in snp_ids:
    record = Entrez.read(Entrez.elink(dbfrom="snp",
                                  id=snp_id.replace('rs',''),
                                  db="gene"))
    results = record[0]['LinkSetDb'][0]['Link']
    for result in results:
        uid = result['Id']
        handle = Entrez.esummary(db="gene", id=uid)
        uid_record = Entrez.read(handle)
        handle.close()
        uid_summary = uid_record["DocumentSummarySet"]['DocumentSummary'][0]
        gene_name = uid_summary['Name']

        print(snp_id, uid, gene_name)

rs8179065 374654 KIF7
rs12899845 374654 KIF7
rs12899845 90381 TICRR
rs8179066 374654 KIF7
rs80357332 672 BRCA1
rs80356920 672 BRCA1

I can show you a simple way in R, using biomaRt. Let's say you have two snp ids you want to exmine.
SNPids <- c("rs431905509", "rs431905511")

library("biomaRt")

# To show which marts are available
listMarts()

# You need the SNP mart
mart <- useMart("ENSEMBL_MART_SNP")

# Find homo sapiens
listDatasets(mart)

# This will be the dataset we want to use
dataset <- useDataset("hsapiens_snp", mart=mart)

# Show available filters
listFilters(dataset)

# Now list all available attributes
listAttributes(dataset)

# To get the ensembl gene id belonging to the SNPs
getBM(attributes=c('refsnp_id', 'ensembl_gene_stable_id'),
      filters = 'snp_filter',
      values = SNPids,
      mart = dataset)

    refsnp_id ensembl_gene_stable_id
1 rs431905509        ENSG00000070371
2 rs431905509        ENSG00000100075
3 rs431905509        ENSG00000260924
4 rs431905511        ENSG00000145335

If you want a fast (binary) search with a locally-stored dataset, where you know exactly what SNPs and gene annotations are being queried, here's a way to set up files to do your own querying, which should scale up nicely if you have a lot of SNPs.
1) Get annotations of interest and write them into a sorted BED file.
For example, here's a way to get Gencode v26 gene annotations:
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_26/gencode.v26.basic.annotation.gff3.gz \
    | gunzip -c \
    | convert2bed --input=gff - \
    > gencode.v26.bed

2) Get SNPs and write them into a text file sorted by SNP ID.
For hg38, for instance:
$ LC_ALL=C
$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/snp147.txt.gz \
   | gunzip -c \
   | awk -v OFS="\t" '{ print $5,$2,$3,($3+1) }' \
   | sort -k1,1 \
   > hg38.snp147.sortedByName.txt

3) Install pts-line-bisect:
$ git clone https://github.com/pts/pts-line-bisect.git
$ cd pts-line-bisect && make && cd ..

4) Run a query:
$ rs_of_interest=rs2814778
$ ./pts-line-bisect/pts_lbsearch -p hg38.snp147.sortedByName.txt ${rs_of_interest} \
   | head -1 \
   | cut -f2- \
   | bedmap --echo --echo-map-id-uniq --range 1000 - gencode.v26.bed \
   > answer.bed

The file answer.bed contains the genomic interval of the SNP of interest, along with the IDs of any Gencode v26 gene annotations within 1000 nt upstream and downstream of the SNP.
To repeat step 4 within a Python script, you can use subprocess.check_output():
...
pts_lbsearch = os.path.join(project_dir, "pts-line-bisect", "pts_lbsearch")
if not os.path.exists(pts_lbsearch):
  print "Content-type:application/json\r\n\r\n"
  print json.dumps({'error' : 'could not locate pts_lbsearch [%s]' % (pts_lbsearch)})
  sys.exit(errno.ENOENT)

if term.startswith('rs'):
  query_fn = os.path.join(annotations_dir, "hg38.snp147.sortedByName.txt")
  cmd = "%s -p %s %s" % (pts_lbsearch, query_fn, term)
  try:
    res = subprocess.check_output(cmd, shell=True)
    if res and len(res) > 0:
      (chr, start, stop) = res.strip().split('\n')[0].split('\t')
    else:
      pass
  except subprocess.CalledProcessError as err:
    sys.stderr.write("Query failed!\n")
    sys.exit(errno.EIO)

#
# map 'chr', 'start', and 'stop' variables to annotations with
# another subprocess call (not shown -- use example above)
#
...

Have you tried using SnpEff and SnpShift (http://snpeff.sourceforge.net) for annotation of the SNPs?
Even if your desired organism available in their annotation database you can create them easily with the gene bank files. Other than annotation the tools features a bunch of other amazing functionalities.

using ncbi-efetch and XSLT with the following stylesheet:
<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet
    xmlns:xsl='http://www.w3.org/1999/XSL/Transform'
    xmlns:r="https://www.ncbi.nlm.nih.gov/SNP/docsum"
    version='1.1'
    >
<xsl:output method="text"  encoding="UTF-8"/>

<xsl:template match="/">
<xsl:apply-templates select="r:ExchangeSet"/>
</xsl:template>
<xsl:template match="r:ExchangeSet">
<xsl:apply-templates select="r:Rs"/>
</xsl:template>
<xsl:template match="r:Rs">
<xsl:apply-templates select="r:Assembly/r:Component/r:MapLoc/r:FxnSet"/>
</xsl:template>

<xsl:template match="r:FxnSet">
<xsl:text>rs</xsl:text>
<xsl:value-of select="../../../../@rsId"/>
<xsl:text>  ncbi:</xsl:text>
<xsl:value-of select="@geneId"/>
<xsl:text>  </xsl:text>
<xsl:value-of select="@symbol"/>
<xsl:text>
</xsl:text>
</xsl:template>

</xsl:stylesheet>

usage (group the rs# by 3):
$ echo -e "rs8179065\nrs12899845\nrs8179066\nrs80357332\nrs80356920\n" | xargs  -L 3 echo  | tr " " "," | while read R; do curl -sL "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=snp&id=${R}&retmode=xml" | xsltproc --novalid  transform.xsl -; done | uniq
rs8179065   ncbi:374654 KIF7
rs12899845  ncbi:90381  TICRR
rs8179066   ncbi:374654 KIF7
rs80357332  ncbi:672    BRCA1
rs80356920  ncbi:672    BRCA1

BiocManager::install('biomaRt')
require(biomaRt)
ensembl <- useMart("ENSEMBL_MART_SNP", dataset = "hsapiens_snp")
getBM(attributes=c("refsnp_id", "chr_name", "chrom_start", "chrom_end",
                   "allele", "mapweight", "validated", "allele_1", "minor_allele",
                   "minor_allele_freq", "minor_allele_count", "clinical_significance",
                   "synonym_name", "ensembl_gene_stable_id"),
      filters="snp_filter", values=c("rs6025","rs431905509"),
      mart=ensembl, uniqueRows=TRUE)


text:
What is the significance of the effect size on WOMAC scores for TKA?
https://www.medicaljournals.se/jrm/content/html/10.2340/16501977-2210
I understand the difference between statistical significance and effect size, but I am having trouble interpreting the data on Table III in the article linked above. I am evaluating the WOMAC scores. My first question is can effect size confidence intervals cross zero and still be valid? My second question is about the negative and positive values of the Hedges'g. The first study scores the WOMAC in a way where the higher the score, the better the outcome. The second two studies score the WOMAC in a way where the lower the score, the better the outcome. Does this mean that in the first study a positive Hedges'g indicate larger effect for the treatment group and a negative Hedges' g for the other two studies indicate larger effect for the treatment group?

Hi and welcome to the site. Please reproduce the relevant part of the paper here so we can understand what you are asking.

higher score = higher function Hedges'g g of -0.63 CI95% (-23.9, 5.9)

lower score = higher function. Hedges' g of 0.18 & CI95% (-2.9, 4.3)

Hedge's g measures the difference between two groups here experimental and  control in terms of combined (weighted) standard deviation. A higher g implies greater difference and lower g means lesser difference between  two groups.
It depends  on the study how the scoring of questions is proceeded with : the effect-size will be same. For example, answer to question item may be scored 1 to 10 depending on whether question item is written in positive  sense (the apple is good for health) or negative sense (the apple is not good for health). Positive Hedge's g does not indicate larger effect for the treatment group. THE negative Hedges' g for the other two studies indicates positive effect for the treatment group if the difference has been computed by  deducting experimental group mean from the control group mean (which is less than exp. Group Mean). Overall, I find that a knowledge of questionnaire scaling and scoring of scales and question items could be very helpful in general and in doing a meta-analysis effectively.


text:
Voom function from limma package and Normalization on counts data
I know that Voom function from limma package from Bioconductor converts raw counts into log-CPM values and then Normalization is applied on that, with normalize.method argument.
I would like to know clearly for RNA-Seq data:

how this normalization is done (statistically) ?

normalization across samples and normalization across genes. But don't the know how it is done statistically.

As you answered yourself @Death Metal, voom will by default not perform additional normalization.
However in virtually all cases you would want to do some kind of normalization at least to correct for differences in sequenced reads between the samples. This is why in the manual (page 71 right at the top) the calcNormFactors function from  edgeR is used. This does a global TMM normalization for differences in library size.
Basically voom assumes that you did a normalization and will not do a second one because this could cause all kinds of problems.
In some cases you may prefer normalization other than TMM, for example if the number of genes detected in your samples is widely different. In this case is could be better to do quantile normalization that is available through the normalizeBetweenArrays function in voom.
However, this may over-correct your data and get rid of many interesting differences. In these cases I have found qsmooth quite helpful. Here quantile normalization is only applied across groups if the local distributions are similar enough.