What is a read count? A very basic question: What exactly is a read count in shotgun sequencing and how is it calculated? I'm particularly interested in understanding the definition of the read count in the context of shotgun metagenomic sequencing though I assume the same concept applies to other types of shotgun sequencing as well. Surprisingly, a Google search did not results in any answers!

In the context of metagenomics, the taxonomic profiling of metagenomic data provides a list of microbial taxa and their "counts". These counts are often normalized to account for genome size or are normalized to convert counts to proportions. Understanding the concept of a read count may help me better understand all these.
Thanks!

To add to the answer from Bioathlete, "reads" in the traditional genomic sense are the string of nucleotide data inferred from the sequencer from a given fragment of genetic material (could be Illumina, Nanopore, PacBio, 10X, etc.).

So read count first applies for all of the "reads" defined above generated by the sequencer. From a NovaSeq, the total number of bases from a single flowcell is up to 3 trillion bases at 150bp each read for the highest throughput flowcell (so a read count of 3,000,000,000,000/150 = 20B total read count).

The read count for a single Nanopore flowcell is roughly (30 billion bases/avg read length), but the read lengths can be SUPER long (100kb compared to 150bp is substantial). Long read technologies are generally used for different analyses than short reads.

So all of the above reads would be defined as "raw reads", meaning they are reads not yet analysed in any way to be used for data analysis. So using illumina as an example, for a single Whole Genome sequencing sample we might have 800,000,000 raw reads. We then use an aligner like bwa mem to map the reads to a single reference genome (or in your case you would use a different aligner to map to multiple genomes). From here we still have the concept of raw reads, but we also now have the total number of mapped reads. A typical good illumina library for a single reference genome will have >90% mapping rate for a human sample.

We can then say we have 800,000,000 raw reads, and 720,000,000 mapped reads with a mapping to a reference genome of 90%. So we can compare "raw read count" to "mapped read count" to get an idea of the quality of our library.

In the case of metagenomics, you would consider your "raw read count" to still be the total number of reads you got off the sequencer, but the "mapped read count" would be the sum of the reads mapped to each individual reference genome to get a total quantity of "mapped reads" which could give you an indication of library quality. I'm not sure what a good metagenomics library "mapped read count" is, but I'm sure you can find that information.

Lastly, we want to get an idea of the depth of coverage for each mapped read to the genome. This will tell us how well our mapped reads actually covered the reference genome or genomes we are interested in.

I hope this was helpful!


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?

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.


RSeQC provided bedfile does not work with read_distribution.py I frequently make use of RSeQC, I used to use it with my own bed file, generated using gtf2bed, but now I thought I'd use the one provided by RSeQC themselves from here, I got the file hg38_Gencode_V28.bed.gz, when looking at it with IGV it looks similar to my self generated bedfile however, read_distribution.py (a script in RSeQC), does not report any reads in any of the categories (CDS_exons, Introns, etc).
Any idea why?
Edit: I see that my own bed file (from the GTF) had ENST identifiers like this ENST00000466557 whereas the file from RSeQC supplies them like ENST00000611746.1 (with decimals) that seems to be the only difference as far as I can see.

Ok, I figured it out already, it's and old issue I come accross every now and then. The original bed file lists the chromosomes as chr1,chr2,...,chrX whereas my bam files uses 1,2,...,X

I fixed the bed file like this (python3):

import pandas as pd
bedFile = '/hg38_Gencode_V28.bed'
bedFile = pd.read_csv(bedFile, sep='\t', header=None)
bedFile[0] = [value.strip('chr') for value in bedFile[0]]
newBedFile = 'hg38_Gencode_V28_reformatted.bed'
bedFile.to_csv(newBedFile, sep='\t', header=None, index=None)

And now read_distribution.py works.


cellranger mkfastq with full path to --id flag I have always used cellranger mkfastq to demultiplex 10x genomics runs manually, though recently the commands to do so have been incorporated into a script that should call the cellranger whenever it detects that a run is finished (and can be demultiplexed).

The main difference when calling manually compared to being called automatically by the script is from where the commands are being called.

So calling the commands manually from the directory where the samplesheet is and where the outputs shall be saved with:

/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=demux_10x --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc

does the job.

but now cellranger is being called (automatically) from /home/ngs. The error message that I am getting from cellranger is:

/home/ngs/installed/cellranger-3.1.0/cellranger mkfastq --id=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ --run=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX --samplesheet=/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/20200507_NovaSeq_129_Chromium_Anne_Fred_Samplesheet.csv --qc
/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/bin
cellranger mkfastq (3.1.0)
Copyright (c) 2019 10x Genomics, Inc.  All rights reserved.
-------------------------------------------------------------------------------

/home/ngs/installed/cellranger-3.1.0/cellranger-cs/3.1.0/tenkit/bin/_mkfastq: line 531: __/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/.mro: No such file or directory
Martian Runtime - '3.1.0-v3.2.3'

Invalid name: /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/ (only numbers, letters, dash, and underscore allowed)

According to the error message I cannot use slashes / in the --id flag.

Does anyone what I should pass to --id in order to get the outputs saved as
/nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/demux_10x/?

I already know that if I just pass --id=demux_10x the outputs will be saved at /home/ngs/demux_10x. (again, I am calling the script that calls cellranger from /home/ngs/)

Obviously I could change to the desired output directory i.e. /nexusb/Novaseq/200507_A00154_0492_AHL72WDRXX/ and do --id=demux_10, though I am trying to avoid that for now.

You should be able to use the --output-dir argument to specify a folder of your choice and use --id=demux_10x as before.


What does genotyping data and beta values for SNPs mean for GWAS study? I have a list of Genes with their alleles (RS IDs) associated with a specific condition. I have to find their genotyping data and the beta values from published studies. I looked through several studies and none of them had anything mentioned about the said terms?
Apparently its for calculating the polygenic score of the listed SNPs. Does that mean the genotyping data is the allelic possibilities and beta score as effect size? Effect size is the population size (n)?

The beta statistic is a common output field produced by programs that carry out genome-wide associations. There's a good explanation here.
I interpret a beta statistic as the degree to which a change in genotype from the reference allele alters the phenotype of interest. The association (and beta statistic) can be negative or positive. It should always be tied to a particular variant, usually called the "risk allele".
The beta statistic is related to the number of individuals found in each group with a particular genotype. In that sense it could be argued to be related to the population size (at least, assuming pre-grouping of populations by genotype), but that's not a typical interpretation.
The GWAS catalogue is a good place to start off looking for data associated with GWAS:
https://www.ebi.ac.uk/gwas/home
However, I think that's for summary statistics. If someone has asked you for raw genotype data, that is trickier to get, and tends to depend on the condition. Genotype data is personally-identifying information, and there are lots of issues associated with a blanket public release of such information, which is why it's so hard to find.
Without more information on the particular phenotype, I can't give much guidance, but useful places to look for most traits are dbGaP and UK Biobank. Both of those data sources need forms filled out and approved (to demonstrate an ethically-responsible research aim) in order to access genotype-level data.
There are a few public datasets around that include data uploaded by individuals (who have explicitly consented to their genotype data being public), but it tends to be trait-specific.


Economist article on coronavirus I am wondering about an article in the Economist here:
https://www.economist.com/briefing/2020/02/29/covid-19-is-now-in-50-countries-and-things-will-get-worse
There is a graph there

The explanation is as follows
The course of an epidemic is shaped by a variable called the reproductive rate, or R. It represents, in effect, the number of further cases each new case will give rise to. If R is high, the number of newly infected people climbs quickly to a peak before, for want of new people to infect, starting to fall back again (see chart 2). If R is low the curve rises and falls more slowly, never reaching the same heights. With sars-cov-2 now spread around the world, the aim of public-health policy, whether at the city, national or global scale, is to flatten the curve, spreading the infections out over time.

I am having trouble with this. There are 2 curves here. First of all, I am not sure of the different processes that determine them. It says one is without measures. The other is with measures. But it also says one is with high R and the other with low R. Is R a fixed constant for a disease? Alternatively, is there some variability in R and does that degree of variation affect transmission? This is not explained or stated in the article and I haven't heard R described in this way. At best it is not written clearly in my opinion.

In addition, does without measures necessarily imply high R? Is that what they are getting at? Also, it looks to me as though the blue graph goes to zero very quickly, basically when t is equal to the mean of the orange curve. That means that without measures the epidemic will end quickly (because it exhausts all the possible patients).

Finally the areas under the two curves is about the same. This means that the number of patients is about the same with or without measures.

Anyway, can anyone clear up any of this?

R is a reference to R0, otherwise known as the basic reproductive rate, and means the number of new cases from a single patient infection. R0 has to be above 1 for a disease to persist. It is a classic equation which account for the probability of transmission against immunity. The y-axis is the number of infecteds, very poorly presented.

The idea is that R0 will be lower if there is containment (less people to infect because you can't move around), but will last longer, possibly because the number of naive patients (not encountered the virus) will be larger over a longer period of time.

That particular graph appears to be saying that everyone at risk of infection will get infected with or without control measures (the surface area under the graphs is approximately the same*) and the only thing that will halt the virus is acquired immunity. If this is mainstream economic thinking then I beg to differ, albeit I note the source is CDC (again see note*).

For example, all other respiratory viral infections are seasonal. If COVID-19  is seasonal containment is the right choice, because the window of transmission is limited.

*, Given this is of CDC origin I think either its a mistake, or at the very least was taken out of context, because the idea is that overall there is a lower number of people infected as a result of "social distancing" and that does not appear to be what the two curves show.


Do you know a program/script to count how many sequences have mismatches on each position of my primers? I'm working with primer analysis and design, and so far I've performing it in Excel but: a) I don't like Excel, b) I do it manually, c) I don't like Excel, d) it gets tedious when working with thousands of sequences and tens of primers. My analysis consists of establishing how many sequences have mismatches on each position of my primers, and also count how many mismatches each sequence has with the corresponding primers (see attached image). This last analysis I've been able to (partially) solve it by using ncbi-blast+ tool. However, I will find very helpful a program/script where you input your sequences and primers fasta files and retrieve the number of mismatches at every primer position. Does anyone know something that could help me achieve all this? If it also has an output tab-separated file as the one attached that would be perfect for easy visualization of results. Suggestion/orientations are also grateful! I just want to know if there's something to help me out, otherwise I'll try writing a script to automatize my analysis. Thank you!

I'm amazed at your tenacity in getting Excel to do what you want; that looks like many hours of work to count mismatches.
The hammer I use at the moment for this sort of stuff is LAST, which will do local alignments of sequences to references. If the number of mismatches are not too large (as would be expected for primer sequences), a local alignment will extend to the end of one of the sequences, and I expect that should work fine for your purposes.
I'm going to demonstrate two options, one that can be done just using Excel and web services, and one that needs command-line access. Note that these are only two ways to fetch mismatch counts; you might want to do more than that (it's not clear from your question), and there are many other ways. Also, this may not work for your particular circumstances, but I think it should be quick enough that it's worth a try:
Option 1 - Web-based LAST
In this, I use the Web version of LAST to generate alignments for the primers [as "reference"] mapped to the sequences [as "query"].

To help create the mismatch (or in this case more correctly, "edit distance") statistic, in the Advanced Parameters section, set the a: gap existence cost to 1, and the b: gap extension cost to 1.
Click on Align!
Once the Alignment Results appear, click on the Alignments in tabular format link down the bottom, and save the file to your computer as last.tab.csv (the .csv at the end is the most important bit, to make sure that Excel recognises it as something it can open).
Open that file in Excel, set the delimiter to <tab>
The first column should be the alignment score, which in this case is the match score. This can be converted into a mismatch score by subtracting it from the primer sequence length (column 6).

Option 2 - Local LAST Installation
Installation
Install LAST and samtools.
Index Generation
The first step is to generate a reference index. LAST doesn't care too much about what's query and what's reference, but because I want all the genome sequences to appear on each line, the genome sequences will be the query, and the primer sequences will be the target. I use Martin Frith's lastdb mantra to do this:
$ lastdb -uNEAR -R01 primers.fa primers.fa

Mapping
The next step is to do the actual alignment with LAST. I set the match and mismatch to be both 1 (i.e. the default for the Web version), set the gap existence and extension to be 1, change the output format to tab, and save the result to last.tab.csv:
$ lastal -r 1 -q 1 -a 1 -b 1 primers.fa sequences.fa -f tab > last.tab.csv

Aggregating data
Finally, it's a matter of aggregating the count information. I use grep -v to exclude lines that begin with a #, use awk to subtract field 1 from field 6, then sort the results and use uniq -c to count the mismatches:
$ grep -v '^#' last.tab.csv | awk '{print ($6 - $1)}' | sort | uniq -c
2118 0
   5 1
   5 2


Alignment with arbitrary number of mismatches or gaps I have 23bp long reads and want to find all possible alignments of them to the human genome (hg19, hg38) for an arbitrary number of mismatches (<7), possibly also small indels. I've read in literature that people use bowtie2 for this, so I've tried the following:

echo "GGAGAACAAGTGGCAGATAGAGG" | bowtie2 --all -f -U- -x hg19 -r --no-hd --quiet --score-min C,-36 -p 3 -L 4

I'm not very experienced with bowtie2, and I hoped that this would give me all alignments of the piped-in sequence with up to 6 mismatches (6*-6 = -36) to the hg19 genome which I previously indexed from fasta files. It however only gives me the correct full alignment:

0   0   chr2    173327389   255 23M *   0   0   GGAGAACAAGTGGCAGATAGAGG IIIIIIIIIIIIIIIIIIIIIII AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:23 YT:Z:UU

I've read that bowtie2 by design only allows a maximum of one mismatch per seed. Even if I implement the smallest possible seed length and interval (-L 4 -i C,1 -N 1), I only get this one result which I think is highly unlikely, given that I'm aligning to the whole genome.

What would be the right way of achieving this using bowtie, or a different tool?

This is not something mappers will be good at, because you're going to hit most of the genome. Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:

for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp


Solutions for managing data in a small bioinformatics / 'omics lab? A different sort of problem: even a small 'omics lab generates a lot of data, raw, intermediate and processed. What (software) solutions exist for managing this data, such that "old" data can be retrieved and checked or re-analysed, even after people have left the lab? Important points would be:

ease of installation
ease of putting data in, in an appropriately tagged / labeled fashion (it's no good if the repository is just a centralised bad mess)
useful search and exploration
security (i.e. restricted to members of the lab)
previews / summaries of the data
can accommodate any dataset
local, not SaaS

It seems from my explorations of the topic that there's very little between the primitive (e.g. handrolled Access databases) and big industrial solutions. Things I have looked at include:

Dataverse: very popular, installation seems complex and unclear if uploading is that easy
DSpace: mostly for publications and documents
CKAN
OSF: used this for a while, integrates with a lot of services but uploading data seemed awkward

I don't know of any prebuilt products, but I can describe how we managed this in my postdoc lab, and how I plan to manage it in my newly-started group.

Rule 1: All work happens in the projects project directory on the central
filestore, not on your desktop or laptop.

Rule 2: Heavy computational work is done by the groups standard analysis pipelines. Interpretation is done in jupyter notebooks or Rmarkdown.

A project has a directory on the group's filestore. That directory has a fixed structure:

proj001----raw_data             *
       |
        --external_datasets
       |
        --src                  *
       |
        --notebooks            *
       |
        --pipeline1
       |
        --pipeline2
       |
        -- etc
 * link to a seperate, backed up filesystem.

Pipelines are where heavy duty analysis happens and everyone uses the same standard pipelines, producing the same standard output files and database structure.

So common pipelines might be the mapping pipeline, the readqc pipeline, the differential expression pipeline, the exome pipeline, the chip-seq pipeline etc...

Pipelines have three important outputs: an automatically generated html report, an sqlite database and files in the export directory.

When we archive a project, this is what we save, along with the pipeline's configuration file and log file.

So if I know that Jane Bloggs did an RNAseq 5 years ago on a cell type i'm interested in, if I konw that that was project 5, I know that in the project 5 directory there will be a diff_expression_pipeline directory, and that it will contain an sqlite database called csvdb and that that will have a table called refcoding_deseq_gene_diff and that table will follow a known format. There will be bigwigs in the export directory. Or the BAMs will be in the mapping_pipeline directory.

Of course the problem remains knowing that Jane Bloggs did this RNAseq and what the project was called. We use a Wiki for this, but its not ideal.


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?

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.


How can I run the latest version of samtools on Slurm? I am trying to run the following script on slurm to extract information from .sam files. While everything is perfect on my local machine with the samtools version 1.11, it does not work on slurm. I have this error:
samtools: command not found

That means slurm cannot understand the module if I am not mistaken. What would be a solution here?
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
# module load samtools/1.11 # this is what I tried first
module spider samtools/1.11 # this does not help either...

for FILE in *.sam;

do

        NUM=`echo $FILE | sed -r 's/.sam//' | sed -r 's/nastya//'`

#         echo "$FILE" ;

        samtools flagstat -O tsv $FILE > file${NUM}.txt

done

I would start with:
module avail

Then look to see (a) if there is a samtools module avilable and (b) exactly what it's called.
If it's called samtools/1.11 then:
module load samtools/1.11

Should work, if it doesn't speak to your sys admin.


Calculate genome coverage and depth from alignment I have a .bam alignment file and a genome reference .fasta file. I am looking for a easy to use tool (that I can reference in a publication) to calculate the percentage coverage of the reference by the alignment and if possible also: average depth of reads that cover the gene, min/max depths. Running plot would be nice but not essential.
I wrote some python code to calculate coverage percent, but I want to reference the tool I use in a paper and dont want to mess with making a git repo and linking to it in the paper just for this simple analysis.

You could use samtools coverage as explained in the manual of samtoools.
Here is a example which is also described on the manual site.
samtools coverage -r chr1:1M-12M input.bam

#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8


SNP vs Point Mutation What is the difference between a Single Nucleotide Polymorphism (SNP) and a point mutation? I am quite confused in understanding these term as both of them refer to one base difference from the reference sequence?

The difference between the two depends on to whom you talk ;)

You are right: both refer to one base difference from the sequence. Usually the term "mutation" is used if the change has an impact on the phenotype. The "P" in "SNP" means "polymorphism". A lot of people use this term if more than 1% of the individual in a given population have this variant.

My recommendation: only use "mutation" if you are sure that your variant has an important impact on the phenotype. Only use "SNP" if you are sure that it is a commonly occurring variant. In any other case just use "variant". It is slowly becoming more common to talk about "SNV" (Single Nucleotide Variant) than about "SNP".


log2fold change as heatmap From a paper:

  (D) Expression analysis of multiple lineage-specific differentiation markers in WT and PUS7-KO EBs (14 days). Heatmap shows log2 fold change (FC) PUS7-KO to WT for each individual gene (rows) in three independent experiments (columns).

They have analysed the data in EdgeR but I was wondering how did they plot fold change when they have two condition but they have three independent experiment, I haven't used EdgeR so far, use DESeq2 where I group my biological replicate. For example if I have two condition with three biological replicate then I would group them and analyse, in case of DESeq2 I would get final result with one fold change report.

How did they do in this case reporting three fold change ?

Data sets where I see three wild type and three knockout.

If I read quickly thru the paper, I see you refer to Fig 5D. In the methods they describe that they use Fluidigm (i.e. RT-qPCR) for gene expression, not RNA-seq. The data link you give is for small-RNA-seq.


Database of copy number alterations in different cancer types Is there a database of somatic copy number alterations (CNAs) in different cancer types? That is, annotations of which CNAs are more commonly seen in different cancer types.

I don't know if there's a database that does exactly what you want, but there are some places that might help you figure this out, especially if you already have a list of CNAs/genes/regions in mind.

ICGC is the first that comes to mind, as it has samples from many types of cancer and copy number data for most of them. Depending on what you want though, you might have to apply for access to download the data you want. It's rather annoying, and their data download tools are atrocious to use, but they have a lot of data.

DGV might be another option, though their data browser is rather dated feeling. They publish a lot of gold-standard germline CNVs and such.

COSMIC is a large database for cancer-associated mutations, including CNVs, so it might be worth a look for you.


Bonferroni count = 1 in GO-Slim overrepresentation test (PANTHER) is it possible? I encountered a problem while performing analysis using the PANTHER website. I run the statistical overrepresentation test for a genelist of 99 genes (GO-Slim BP, Binomial, Bonferroni) in February 2020. My result was Bonferroni count = 1 and a list of several enriched GO terms. Next we run the same analysis with the same settings in July 2020 (two days ago) and the result is dramatically different, the Bonferroni count is 1770 and no enriched terms are found.
I am confused - is the result Bonferroni count = 1 possible? Why these results are so different?
This is how the correction is described on PANTHER website:

The expression data analysis statistics now include a Bonferroni correction for multiple testing. The Bonferroni correction is important because we are performing many statistical tests (one for each pathway, or each ontology term) at the same time. This correction multiplies the single-test P-value by the number of independent tests to obtain an expected error rate.
For pathways, we now correct the reported P-values by multiplying by the number of associated pathways with two or more genes. Some proteins participate in multiple pathways, so the tests are not completely independent of each other and the Bonferroni correction is conservative. For ontology terms, the simple Bonferroni correction becomes extremely conservative because parent (more general) and child (more specific) terms are not independent at all: any gene or protein associated with a child term is also associated with the parent (and grandparent, etc.) terms as well.
To estimate the number of independent tests for an ontology, we count the number of classes with at least two genes in the reference list that are annotated directly to that class (i.e. not indirectly via an annotation to a more specific subclass).

GO-Slim BP, Binomial, Bonferroni

This seems easy you because I suspect have run too many tests at 5%.
If you run 10 tests at 5% rejection threshold, the Bonferroni correction threshold is 0.5% - you divide the number of tests by the rejection threshold.
One possibility is the rejection threshold is given in number of genems: in one instance you've run a pairwise test with several thousand genes so the rejection threshold has been set accordingly. The other test you ran I suspect was a multi-variate test, thus the shift in rejection threshold measured in number of genes.

The altnative is that Bonf = 1 is in fact P = 1, in other words there is no significant difference at all.  If P = 1 then this could be no genes show deviation hence 1770, because 1770/1770 =1
Looking at it again the above seems a better explanation.


Quantifying reads mapping to multiple loci I have been using STAR for our RNA-Seq samples. The final.out log file reports percentage of uniquely mapped reads along with percentage of reads that map to multiple loci (less than or equal to 10) and percentage of reads mapping to too many loci (greater than 10). However, I want to break down the multiple loci part to individual counts: Reads mapping to 2 locations, 3 locations, 4 locations .. 10 locations.

The NH tag seems to be used by STAR. However a naive read counting approach results in it reporting more number of reads than total reads.

For example, my final.out looks like this:

   Mapping speed, Million of reads per hour |       1403.36

                      Number of input reads |       53015978
                  Average input read length |       26
                                UNIQUE READS:
               Uniquely mapped reads number |       368916
                    Uniquely mapped reads % |       0.70%
                      Average mapped length |       26.45
                   Number of splices: Total |       1057
        Number of splices: Annotated (sjdb) |       0
                   Number of splices: GT/AG |       802
                   Number of splices: GC/AG |       1
                   Number of splices: AT/AC |       0
           Number of splices: Non-canonical |       254
                  Mismatch rate per base, % |       0.31%
                     Deletion rate per base |       0.00%
                    Deletion average length |       1.45
                    Insertion rate per base |       0.00%
                   Insertion average length |       1.00
                         MULTI-MAPPING READS:
    Number of reads mapped to multiple loci |       45766732
         % of reads mapped to multiple loci |       86.33%
    Number of reads mapped to too many loci |       3757890
         % of reads mapped to too many loci |       7.09%
                              UNMAPPED READS:
   % of reads unmapped: too many mismatches |       0.00%
             % of reads unmapped: too short |       5.89%
                 % of reads unmapped: other |       0.00%

Counting histogram of number of positions a read maps to using pysam:

def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts

results in

Counter({1: 330606,
     2: 86772164,
     3: 329,
     4: 38083,
     5: 31,
     6: 1094,
     7: 129,
     8: 50,
     10: 50})

The count 1 reads are fine even though they do not match the counts in final.out file since I am counting a certain category of reads (say those mapping to tRNA only), but the reads mapping to 2 locations are highly overestimated. Why is that?

You almost had the correct python code already, you just need to filter out secondary alignments:

def get_reads_hist(bam):
    bam = pysam.AlignmentFile(bam, 'rb')
    counts = Counter()
    for query in bam.fetch():
        if query.is_secondary:
            continue
        nh_count = Counter(dict(query.get_tags())['NH'])
        counts += nh_count
    return counts


Explanation for RNA-seq samples not clustering in PCA as expected A colleague is analysing RNA-seq data - the study design is 2 treatments, 3 replicates, 3 tissues. In their PCA plot the samples clustered neatly by tissue. Except for two samples - two tissue samples originating from the same animal clustered with the wrong tissue.
I suggested that the two samples had likely been switched accidently somewhere along the way. They think its explained simply by biological variation which I'm not convinced of. Are there any technical explanations for why samples could cluster on a PCA in such a way?

It's hard to say what is going on without looking at the plots. But it is a reasonable guess that samples are mislabelled and I think it is worth checking the experimental log carefully. It is also possible that the two samples are of low quality. People often use PCA to check and remove outliers from their RNA-seq libraries. So you may also check the RNA-seq library quality for those two samples.


How to visualize genome track of gene in specific cell-lines? I'm trying to make a plot showing genome tracks of specific genes in specific cell-lines of RNA-seq and Chip-seq data.

It should look something like this I have recently seen this encode, but in this I didn't find any MDAMB231 cell-lines.

May I know how to make such plots.

I got the Counts data from CCLE portal.

             BT-549 MDAMB231
RP11-554I8.2    0   246
LINC00707      50   107
   RBM17    17036   10819
  PFKFB3    18098   7668
RP11-393K12.2   0   0

You can use screen.encodeproject.org
It collected a lot of data from multiple cell lines or tissues. And you can export tracks to UCSC browser or other browsers.


Understanding DESeq2 design, contrast and results I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.
The read counts for the genes are summarized in a file that I load as follows in R:
> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
 [1] "WT_RT_1"          "WT_HS30_1"        "WT_HS30RT120_1"   "prg1_RT_1"
 [5] "prg1_HS30_1"      "prg1_HS30RT120_1" "WT_RT_2"          "WT_HS30_2"
 [9] "WT_HS30RT120_2"   "prg1_RT_2"        "prg1_HS30_2"      "prg1_HS30RT120_2"

I describe the experiments as follows:
> col_data <- DataFrame(
    geno = c(rep("WT", times=3), rep("prg1", times=3), rep("WT", times=3), rep("prg1", times=3)),
    treat = rep(c("RT", "HS30", "HS30RT120"), times=4),
    rep = c(rep("1", times=6), rep("2", times=6)),
    row.names = colnames(counts_data))
> col_data
DataFrame with 12 rows and 3 columns
                        geno       treat         rep
                 <character> <character> <character>
WT_RT_1                   WT          RT           1
WT_HS30_1                 WT        HS30           1
WT_HS30RT120_1            WT   HS30RT120           1
prg1_RT_1               prg1          RT           1
prg1_HS30_1             prg1        HS30           1
...                      ...         ...         ...
WT_HS30_2                 WT        HS30           2
WT_HS30RT120_2            WT   HS30RT120           2
prg1_RT_2               prg1          RT           2
prg1_HS30_2             prg1        HS30           2
prg1_HS30RT120_2        prg1   HS30RT120           2

I want to build a DESeq2 object that I could use to either:

find differentially expressed genes when the treatment varies for a given fixed genotype

or:

find differentially expressed genes when the genotype varies for a given fixed treatment

In the bioconductor help forum I think I've found a somewhat similar situation, and I read the following:

Try a design of ~ genotype + genotype:condition
Then you will have a condition effect for each level of genotype, including the reference level.
You can constrast pairs of them using the list style of the 'contrast' argument.

However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.
So I tried the following more symmetric formula:
> dds <- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)

Now, can I for instance get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"?
And how?
If I understand correctly, the above-mentioned "list style" uses names given by the resultsNames function. In my case, I have the following:
> resultsNames(dds)
[1] "Intercept"               "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120"   "genoWT.treatRT"

I guess I would need a contrast between "genoprg1.treatRT" and a "genoprg1.treatHS30", but these are not in the above results names.
I'm lost.

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:

dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)

The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").

As an aside, this is probably a case where the edgeR-preferred way of creating groups of genotype-treatment combinations and then using a model of ~0+group might make your life a bit easier. You'll get the same results (more or less) regardless, but it'll probably be easier for you to think in those terms rather than remembering that the base level will be treatment HS30 and geno prg1.


How to use extractfeat? I am working on a script on bash linux to get the CDS (coding sequence) of a gene using extractfeat by EMBOSS. It gives me the error: Warning: No sequences written to output file

I am unsure if my syntax is right, but this is what I have:

extractfeat inputFile outputFile -type CDS

Perhaps,removing introns (if it is a higher eukaryote)

extractfeat tembl:x65921 -type CDS -join

It is difficult to know without knowing the gene.

BTW NCBI's Entrez "efetch" will do essentially the same thing, but its more complicated


Reads mapped to exonic, intronic and intergenic regions After the alignment step I checked the rnaseq metrics of all the samples. Among 40 samples three samples show high percentage of reads mapped to intronic regions. What could be the reason?

Samples   Exonic            Intronic    Intergenic
Sample1 545479 (12.8%)  2512309 (58.8%) 1217201 (28.5%)
Sample2 372836 (8.3%)   2556934 (56.8%) 1573032 (34%)
Sample3 529618 (7.3%)   3934615 (54.5%) 2750720 (38.1%)

I also see that 35-40k genes among 58k genes having "zero" read counts. Is this due to contamination? What could be the reason? What does reads mapping to exonic, intronic, intergenic tell?

Update: I used hisat2 for alignment. I used human genome "grch38_snp_tran" from hisat2 website. Libraries were generated using ribosome depletion kit.

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.


vcftools: indel size histogram command returns empty file I would like to get some summary statistics on a vcf file from one individual, which has over a million variant calls. I've tried to make a histogram of indel sizes with this command,

vcftools --vcf sample.vcf --out sample --hist-indel-len

but it returns an empty file. I see thousands of deletions in the vcf, although strangely no insertions. A typical line in the vcf file I have looks like

NC_025968.3 168363 191 N <DEL> . . SVTYPE=DEL;STRANDS=+-:19;SVLEN=-1026;END=169389;CIPOS=-9,330;CIEND=-152,9;CIPOS95=-3,61;CIEND95=-32,3;IMPRECISE;SU=19;PE=19;SR=0 GT:SU:PE:SR ./.:19:19:0

Update: I aligned reads to my reference with BWA MEM, then used samblaster to extract discordant pairs and split-reads. I sorted and indexed the three bam files, then used lumpyexpress to generate the vcf.

samblaster -i sample.sam --excludeDups --addMateTags --maxSplitCount 2 --minNonOverlap 20 | samtools view -S -b - > sample.bam

samtools view -h sample.bam | samblaster -a -e -d sample.disc.sam -s sample.split.sam -o /dev/null

lumpyexpress -B ../sample.sorted.bam -S ../sample.split.sorted.bam -D ../sample.disc.sorted.bam -o sample.vcf

I would like to see the size distribution of deletions, and if possible, figure out why there are no insertions.

It looks like lumpy looks for breakpoints in your reference. Presumably places where the is a haplotype switch detectable with the data on-hand?

My only thought is that since lumpy is being used to produce the VCF file that it is inherently not looking for insertions. I would suggest running the same BAM file through the GATK best practices pipeline (easier said than done) and vcftools pipeline. If those produce insertions but lumpy does not, it must be due to how lumpy produces the calls. Good luck!


Program to make a haplotype network for a specific gene Currently my lab manually uses pop art to make the network after using PHASE to id an individual's alleles. Example:

I just feel like there's got to be a better, more time effective way to do this.

This is a haplotype map, where each node is proportional to the frequency of a given allele. The formal tool to do this is Network 5 found here. The problem with manaul haplotype maps is the area of the node should be proportional to the frequency, if this is not accurate it could be perceived to be misleading.

The approach for looking at each parsimonious state along a branch is a phylogenetics method which is performed via MacClade 4.0 ... requires an old Mac though.


How to find hairpin loop structure in a large set of long non coding RNA transcripts I have done de novo assembly of pair end raw read sequences and resulted transcripts sequences were separated based on coding potential into two categories: long non-coding RNA transcripts and coding RNA transcripts.

I want to find how many of long non coding can make hairpin loop structure and how can I annotate these miRNA which are generated from long non coding RNA in plants.

For the first question: "how many of long non coding can make hairpin loop structure?"

The easiest thing to would be to run your sequences through an RNA secondary structure prediction tool. There's a few tools for doing this but the most commonly used are RNAfold from the ViennaRNA package and mfold.

This will give you an output that looks like this, where the parentheses indicate which bases are paired with each other.

GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA
(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).

To see if there's hairpins, you'll need to annotate the elements within it. Take a look at this question on biostars for an example of how to do this. In short, you can use the forgi library to produce an annotation like this:

(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).
sssssssmmsssshhhhhhhhhssssmssssshhhhhhhsssssmmmmmssssshhhhhhhsssssssssssse

If there's an h in the output annotation, you have a hairpin. The remainder of your question about annotating miRNAs will be left to somebody else.


determine if ChIP-seq peaks are broad or narrow Is there a method to determine if the peaks are broad or narrow? ENCODE provides some guidelines:

Although those cover common histone marks, there are many others. If you are using one of the ones that are not listed, is there a way to determine the appropriate peak width?

This question was inspired by this previous post.

I'd say mostly this is a question of understanding the underlying biology and the relevant literature.

If it is not known in the literature whether a mark is peaky or broad, evidence might come from FISH studies or Low throughput qPCR. Another way to look would be to examine the signal expressed as fold enrichment over input on a genome browser and look if peaks tended to be small (a few hundred bp, e.g. a transcription factor), or large (several 1kb, e.g. H3K27me3) or contain both (e.g. RNA polymerase II).


Bacterial DNA at the tail of transcriptome reads. What does that mean? I am assembling a transcriptome obtained from the Internet. The transcriptome was extracted from a human cancer tissue that had been previously grafted into a mouse. I have detected that many assembled RNA sequences do not match precisely to their reference but have additional sequences at one of their tails containing bacterial DNA. It's not that there are reads pertaining to that bacterium (always the same one) but that the RNA sequences of the cell have incorporated fragments of bacterial sequences (the reads corroborate these insertions with high coverage), therefore I can't just filter out bacterial reads as there are not reads exclusively pertaining to bacteria but hybrids.
My question is how can it be possible? What possible mechanisms could have brought these bacterial fragments to the RNA of the cell. I have read in some papers that there is speculation that bacteria could be involved in the genesis of cancer cells. That would be an explanation if it wasn't because there are also mouse RNA for which the same circumstance applies: they also have bacterial sequences at the tails. I have also read in papers that there is speculation with the possibility that cancer cells are more susceptible to get DNA from bacteria transferred to their DNA. I am speculating that is the reason in this case (and somehow that weakness was transferred to contiguous mouse cells) but I would like to know if there is another possibility I am missing.
Thanks.

From time to time these situations are encountered, notably when PCR is being deployed in rare DNA samples.
In this instance the wiki states,

Ralstonia solanacearum is an aerobic non-spore-forming,
Gram-negative, plant pathogenic bacterium.".

I'd leave it. The simple explanation is there has been a bad mix up. The best explanation is the people loading the Illumina machine, designating the adaptors etc ... In this scenario it has nothing to do with bacterial genome integration, plant bacteria don't touch mice and aerobic bacteria wouldn't grow in cancer cells.
The second issue is what is a soil bacteria doing in a mouse cancer biological facility? For the staff this wouldn't be good news.
The singular biological  possibility is its a Pseudomonas  and some species are serious opportunistic pathogens of immune-weakened humans, but there are a lot of Pseudomonas species and this one doesn't ring any bells, BUT its mice so who knows? One species in particular is associated with hospital outbreaks, but a biomedical unit that isn't quite the same thing.

The pathogen in humans is called Pseudomonas aeruginosa. If hospitals get this infection it is BAD NEWS, for example burns victims are really susceptible and without intervention it kills them. If you have evidence that the bacterial equivalent in mice is the bacterial pathogen you have found that is a substantial step. Biological Safety Facilities (BSFs) DEFINITELY get infected. They hate it when this happens and can be very picky for fear of a contamination, but it certainly happens. If you have evidence of a major contaminent pathogen in the long run it would be of HUGE help to BSFs globally. No-one really bothers with the disease pathology of mice - because their mice and spend their entire lives going 'squeek' - but when they are models for cancer their status changes beyond recognition. So I'd erase my initial thoughts and suggest, what you have found is publishable and the application is BSF contamination - don't say that though say this could be the mouse equivalent of P. aeruginosa. Its certainly a publication and goodluck.


Recommended workflow for RNA long read Kallisto-like TPM estimation? On short reads, Kallisto/Salmon is a standard workflow for measuring RNA transcript counts in TPM. However, when I tried to Google a similar workflow for long reads, it's not clear there is a replacement.

Stringtie now does long reads too. We've found flair to be quite difficult to install and to crash often. That said, stringtie isn't very good about showing you its working (e.g. I have no idea how to get it to give the assembly it uses to figure out all the isoforms). Its values for TPM/FPKM also seem a little weird, in that they're perfectly correlated with each other, and for transcripts they're perfectly correlated with coverage too.
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1910-1
https://github.com/gpertea/stringtie
TALON is another tool that does isoform assembly. I'm not sure whether it does quantification though.
https://github.com/mortazavilab/TALON#


Retrieve abstract/summary for NCBI bookshelf via Entrez I want to retrieve abstracts/summaries for NCBI bookshelf entries, eg: "NBK1440"

Docs say the dbname is "books" and efetch guide says a rettype of "docsum" works for all databases. However when I make the query:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&rettype=docsum

returns:

<ERROR>UID=1440: cannot get document summary</ERROR>

The default query is useless (just echos back the ID):

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1440&retmode=text

returns

1440

Picking another 3:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=books&id=NBK1391,NBK410087,NBK1391&rettype=docsum

Fails the same way:

<ERROR>UID=1391: cannot get document summary</ERROR>
<ERROR>UID=410087: cannot get document summary</ERROR>
<ERROR>UID=1391: cannot get document summary</ERROR>

The same goes for esummary:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=1440

Dies on NBK1440 and removing NBK gives:

<ERROR>UID=1440: cannot get document summary</ERROR>

The values come from the (9Mb) ClinVar citations file:  https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/var_citations.txt - the first is from tihs line:

15048   9   1800562     NCBIBookShelf   NBK1440

ANSWER

Thanks to arupgsh - yes the ClinVar citation_id was a RID for bookshelf not an ID like it is for Pubmed and PMC. FYI my solution was:

def get_summary_for_bookshelf_rid(bookshelf_rid):
    handle = Entrez.esearch(db="books", term=bookshelf_rid)
    search_results = Entrez.read(handle)
    handle = Entrez.esummary(db="books", id=','.join(search_results['IdList']))
    results = Entrez.read(handle)

    for r in results:
        if r["RID"] == bookshelf_rid:
            return r

    return None

get_summary_for_bookshelf_rid("NBK1440")

returns:

{'PubDate': '2000/04/03 00:00', 'Title': 'HFE-Associated Hereditary Hemochromatosis', 'Text': '', 'RType': 'chapter', u'Id': '1475938', u'Item': [], 'Book': 'gene', 'Parents': '', 'RID': 'NBK1440', 'ID': 'gene/chapter/hemochromatosis/PMC'}

Entrez requires the unique identifiers(UID) for fetching related info. The id you are using in the query is RID. I guess that why you are getting <ERROR>UID=1440: cannot get document summary</ERROR>.

List of example uids using the term "ncbi+blast":

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=books&term=ncbi+blast

Summary of the first item:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=4459608

Using the Entrez python package:

from Bio import Entrez
Entrez.email = "youremail@domain.com"
handle = Entrez.esummary(db="books", id="4459608,4456634", retmode="xml")
records = Entrez.parse(handle)
for record in records:
     print(record['Title'])
handle.close()

I'm not sure about the availability of book abstracts.


For phylogenetic tree construction from core-genome which one is preferable: amino-acid based MSA or nucleotide based MSA? The genomes are from same species.
Is it true that, in phylogenetic tree constructed from amino-acid based MSA (multiple sequence alignment) some information are lost, so for  phylogenetic reconstruction of closely related sequences nucleotide based MSA  should be used?

For a species complex nucleotide phylogenies all the way.

The reason is neutral mutation, which is not observed at the amino acid level because these mutations are part of the degenerate code which don't induce amino acid substitutions. Such mutations are commonly referred to as the third codon wobble, although occasionally the first codon can mutate and not produce a change in amino acid.

Neutral otherwise known as synonymous mutations occur at a much higher rate than mutations that cause a substitution of an amino acid. This is the reason why  you were told information is lost at the amino acid level because you're losing neutral mutations.

There are circumstances when the amino acid phylogenies are preferred to nucleoside trees, but these are bit pretty specialised areas.


RPGC normalisation creates artefacts at centromere I am studying ChIP-Seq data in HeLa cells and I've started using the RPGC normalisation of deepTool's bamCoverage. MACS2 also uses this normalisation for its peak calling.

I am seeing a large number of peaks (like, over 50% of my peaks) at the centromere of each chromosome, and the RPGC normalised BW files also show this increased signal at the centromere.

Does this mean that my protein really is at the centromere of each chromosome, or is it some normalisation bias I'm not aware of?

Here's the general aspect of my peak calling (in green). Notice the large increase at the centromere of each chromosome.

And here's a zoom on the centromere of chromosome 6. In blue at the top is raw BW, then there's a lane for input normalised BW but not RPGC, then green is the peak calling by MACS2, and black is RPGC normalised BW.

UPDATE

Here is some other protein ChIPs that I ran through the same normalisation:

As you can see, they don't seem to be enriched at the centromere. So unless anyone brings other contributions, I will accept this as answer.

RPGC normalisation doesn't appear to create artefacts at centromere. The peaks observed here are as real as can be.


Why is my STAR reference genome indexing aborting on my GNU/Linux server but not on my Mac OS X laptop? I am running the following command to index my genome:

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gff3 --sjdbGTFtagExonParentTranscript Parent --sjdbOverhang 100

and I get the following output:

Feb 27 10:27:37 ..... started STAR run
Feb 27 10:27:37 ... starting to generate Genome files
Feb 27 10:27:40 ... starting to sort Suffix Array. This may take a long time...
Feb 27 10:27:42 ... sorting Suffix Array chunks and saving them to disk...
Feb 27 10:28:26 ... loading chunks from disk, packing SA...
Feb 27 10:28:30 ... finished generating suffix array
Feb 27 10:28:30 ... generating Suffix Array index
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc

Searches for the GNU/Linux error lead to some pages about C++ that I don't understand.

Submitting it as a batch job works fine but I would like to know why it does not work on a node. It is quite practical for me to test small code like this without submitting a batch job.

The command above works fine on my laptop running Mac OSX.

The STAR version is STAR_2.5.3a on both computers: STAR_2.5.3a

The memory usage of the node is the following:

total            used       free     shared    buffers       cached

251              251         0         0       0             232

so, is there 0 memory for me to use on this node?

Why is my indexing aborting?

A std::bad_alloc error tends to mean that you've run out of memory. It's not unusual for head nodes to be fairly limited in capacity, so presumably that's the issue.


Impute phenotype under some constraints Errors happen frequently in the lab and I got one sample that we mixed one important information (the region of the sample) between two samples. Now we don't know which sample is from which region, but we know all the other variables (Date of the sample, patient, and other variables).

So far, my solution has been to make that values as NA to be able to analyze it, but I would like to know if there is any solution. If I can convey to test or other analysis downstream that sample A and sample B are from Location1 or Location2 and if A is from Location1, B is from Location2 and otherwise round.

There are several packages for imputing quantitative data, but I couldn't find any to store this kind of information or to impute qualitative data. I can recall a recent paper with single-cell sequencing that estimated phenotype and some other variables but I can't find it neither I think it is possible with RNAseq.

Although maybe I could build a classifier to actually do that (but that would be too long and I'm not sure if it is worth it)

Notes

My samples are from more than two locations.
My samples consist of 16S microbiome DNA-seq and RNA-seq

Question 1: Is there a tool for imputing phenotypes while also using genetic data?

Yes, try Andy Dahl's Phenix from the Marchini lab.

Disclosure: I did my PhD in the Marchini lab, although I was not directly involved in this research.

Question 2: Is there a way to impute categorical data using methods intended for continuous data?

If the variable of interest has only two categories, then you should be able to simply recode the categories as arbitrary (differing) numbers and get sensible results. But do make sure you understand the tool you are using to see if this will work.

For more than two categories, if your imputation approach allows the simultaneous imputation of multiple variables, you could try recoding your categories as contrasts or one-hot.  After recoding, you could impute the variables resulting from recoding together.


Calculating average coverage for .bam files (sequence data) (Full discolosure that this is my first time working with sequence data, and with the bash scripting.)
I need to calculate the average coverage for any .bam file.
After some searching I wrote the following script:
# Script to calculate the average coverage of a genome sample

# compute the total length of the sample
tot=$(samtools view -H $1 | grep -P '^@SQ' | cut -f 3 -d ':' | awk '{sum+=$1} END {print sum}')
echo $tot

#compute the coverage at each point
samtools depth -ao cov_$1 $1 #your .bam file has depth calculated for each point

#calculate the average coverage
sum=$(awk '{sum+=$3} END {print sum}' cov_$1)
echo $sum
echo

avg=$(echo "$sum/$tot" | bc -l)
echo ''The average coverage of the sample $1 is $avg x.''

The logic behind the calculation of tot is as follows:
samtools view -H  gets the header of the bam file
grep -P gets all lines in this header with SQ
cut -f 3 -d ':' removes fields that are lists of more than 3 letters while -d uses delimiter ':' instead of tab and this leaves only the numbers i.e. lengths of the chromosomes
awk '{sum+=$1} END {tot=sum}' calculates the length/total bases in the chromosome
However, my output seems incorrect. Any advice on where I might be making an error (in code or thought process), or a better method would be appreciated.
One thing that worries me is that my tot variable is the total number of bases, so if the samtools depth function does not calculate coverage for each base then my sum and tot will not be in the same units and my calculated coverage will be wrong...

For a quick estimate you’re making it more complicated than necessary.

The theoretical average coverage is $\frac{n \cdot \hat l}{N}$ where $n$ is the number of reads, $\hat l$ is the average read length and $N$ is the genome size.

samtools idxstats gives you the chromosome lengths and number of mapped reads in one convenient list.

Putting this together, we get
samtools idxstats data.bam \
| awk -vreadlen=120 '
    {
        len += $2
        nreads += $3
    }
    END {
        print nreads * readlen / len
    }
'

(In this case I’ve used 120 as the read length.)
However, note that this over-estimates empirical sequencing depth, mainly due to the fact that a given read may not be fully mapped; thus the the average read length is an upper bound on the average mapped length. If you have an estimate for the average percentage of mapped bases per read, you can multiple that by the above coverage estimate for a more accurate number.

If you need a more accurate estimate, your approach is fundamentally sound (and the result is probably correct!). For low-coverage data you can speed up your analysis by dropping the -a argument to samtools depth: you don’t need to record zero-coverage bases, they don’t contribute to the total. You also don’t need to store the coverage into an intermediate file, thus reducing IO cost:
sum=$(samtools depth "$1" | awk '{sum+=$3} END {print sum}')

Some more remarks on your code:

The -P flag for grep is not supported everywhere and thus limits portability, so you want to avoid using it if possible. It turns out that you don’t need it here: grep '^@SQ' works just as well.
While the philosophy of UNIX is based on solving problems by combining individual tools that each solve one part, it is sometimes convenient to merge individual steps:

You can often combine cut and awk, in this case by setting the field separator:
awk -vFS=: '{sum+=$3} END {print sum}'

You can often combine grep and awk:
awk -vFS=: '/^@SQ/ {sum+=$3} END {print sum}'

Always double-quote your shell variables, otherwise your script will break if the values contain unexpected spaces or other special characters. This is especially important for filenames, because spaces and other special characters are valid characters in paths!
Your code is clean and readable. Well done.


How to get % similarity between strains and mutation files I’m very new to python, and having some difficulty getting hang of some more complicated things
I have multiple files which look like so:

hCoV-19/Singapore/4/2020|EPI_ISL_410535|2020-02-03

hCoV-19/USA/WA13-UW9/2020|EPI_ISL_413601|2020-03-02

hCoV-19/USA/WA-UW142/2020|EPI_ISL_416680|2020-03-11

Please be aware that the lines above are meant to be one file
I want to extract the EPI_ISL_000000 for an easy comparison among files.
Could someone please advise on:

A programme to extract this data into new files (There’s many lines in each file - 1000+)

A programme to then give a % comparison between two or more files - comparing all lines in one file against all lines in a second+ file

left_lineagelist = [x.split('_')[-1].split('|')[0]
                          for x in left_lineagelist]
        right_lineagelist = set([x.split('_')[-1].split('|')[0]
                          for x in right_lineagelist])

Allows for extraction of 6 digit EPI, provided the file has had sequences removed prior; as such:
for line in lines:
    if line[0] == '>':
        print(line[1:])


Using Python, how to convert a pandas DataFrame into a VCF? Let's say I have a pandas dataframe with fields CHROM, POS, ALT, REF. In this special case, I also wouldn't care about IDor FILTER, INFO could be blank (or meaningless) and we'll write QUAL as each 40.

import pandas as pd
example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)
print(df)

   CHROM      POS   REF ALT
0     20    14370     G   A
1     20    17330     T   A
2     20  1110696     A   G
3     20  1230237     T   C
4     20  1234567  GTCT   G

In this case, saving the output with the extension .vcf via df.to_csv("myfile.vcf") with the other required fields would created the tab-delimited format of the VCF, but it wouldn't include the mandatory header. https://samtools.github.io/hts-specs/VCFv4.1.pdf

##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
...
#CHROM POS ID REF ALT QUAL FILTER INFO

The most popular package to work with VCF's appears to be PyVCF: https://pypi.python.org/pypi/PyVCF

However, one needs to begin with a VCF as a template, and add to this (from the documentation):

vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
for record in vcf_reader:
     vcf_writer.write_record(record)

It appears the only information necessary from vcf.Reader(filename='vcf/test/tb.vcf.gz') is the metadata?

print(vcf_reader.metadata)
OrderedDict([('fileformat', 'VCFv4.0'), ('fileDate', '20090805'), ('source', ['myImputationProgramV3.1']), ('reference', '1000GenomesPilot-NCBI36'), ('phasing', ['partial'])])

Otherwise, tb.vcf.gz already has variants, which I think means I would need to feed in a VCF with no variants?:

for record in vcf_reader:
     print(record)

Record(CHROM=20, POS=14370, REF=G, ALT=[A])
Record(CHROM=20, POS=17330, REF=T, ALT=[A])
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Record(CHROM=20, POS=1234567, REF=GTCT, ALT=[G, GTACT])

You will have to define the header either from an existing VCF file or hardcoded into you python script.  Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.

You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ID, FILTER, INFO and QUAL columns

import pandas as pd

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']]

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

Then myfile.vcf looks like

##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40
20  17330       T   A   40
20  1110696     A   G   40
20  1230237     T   C   40
20  1234567     GTCT    G   40


Biopython Failure to break out of loop After qblast, I wanted to print only details from the first 3 BLAST output. Yet, the code below printed the details of every BLAST output. Why can't I break out of the loop to achieve my desired output?

for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
        break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

Assuming that your actual code is this:

for alignment in blast_records.alignments:
    countResults = 0
    for hsp in alignment.hsps:
        if countResults == maxResults:
            countResults = 0
            break
        else:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)
            countResults += 1

Here you have an outer loop, which is iterating over all the possible database matches, and an inner loop, iterating over the alignments within that database entry. You are currently limiting the inner loop, but not the outer loop, which you probably want to do. Try this:

# Take the first 3 database hits
for alignment in blast_records.alignments[:3]:
    for hsp in alignment.hsps:
            print(alignment.title)
            print(alignment.length/ len(seq)*100)
            print(hsp.expect)
            print(hsp.identities)


Get start and end coordinates per chromosome I will use GATK for SNP calling (HaplotypeCaller). I need to feed the interval file in the command, otherwise I get errors (even though I want to use the whole genome, not a subset  - btw it's not the point of the question, so will not elaborate).

The interval file for GATK can be for example a bed file, with 3 columns: chr_name  chr_start  chr_end. I don't have this file, but have the genome and the reads.

To obtain the intervals per chr/scaf, I proceeded this way:

converted the alignment to interval file, using bam2bed    (BEDOPS)
extract the columns $1,2,3 (using cut)
get the min & max coordinates per chr/scaff with a small R script:

bd <- read.table('all.bed', h=FALSE)
bd_min <- aggregate( bd$V3 ~ bd$V1 , bd, function(x) min(x))
bd_max <- aggregate( bd$V2 ~ bd$V1 , bd, function(x) max(x))
bd_coord <- merge(bd_min, bd_max, by='bd$V1')
write.table(bd_coord, file='coords.bed', sep='\t', quote=FALSE, col.names=FALSE, row.names=FALSE)

These should thus represent the mapped intervals, conceptually same as not feeding the intervals at all (= considering all the sequences).

The ideas for this:

in the bed file: $1 is be the feature, $2 the start position, $3 the end position
so, I need to get the min of $2 for each unique item in $1(for the start coord) and the max of $3 for each item in $1 (for the end coord)

The output looks plausible:

$ cat coords.bed | head -3
Bla_chrm1   678 43860826
Bla_chrm10  181 20381540
Bla_chrm11  343 20367560

My question here:

is this a correct way to proceed?
is there a standardized way to perform this?

My main concern is that I will have to re-build this "coordinate" files for every GATK run, because I think the coordinates per chr/scaf will shift, even slightly, at each GATK run with different datasets.

My other minor concern, is that this approach is a bit slow (given the 160M lines in the current bed file), so an unix-tools solution will also be accepted. I tried to compute myself the max/min using awk (e.g. awk '$3>max[$1]{max[$1]=$3; row[$1]=$0} END{for (i in row) print row[i]}'), but I get different results than with the R approach...

Not sure if I completely understood your question, but if you want simply an interval of your entire genome  you could do instead:

samtools faidx MyGenome.fasta

Obviously that means that you need access to that, not sure if you have.
The resulting fai file contains a format similar to that:

000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

Where the first column is your chromosome and the second one your length.
I guess you could then just do a awk '{OFS=FS="\t"}{print $1,"0",$2-1}' to get your file.


What are the largest fully publicly available molecular QTL datasets? I'm familiar with the sQTL and eQTL GTEx data that can be downloaded from the
GTEx Portal, but I'm interested in gathering at least two other types of molecular QTL datasets:

Larger eQTL and/or sQTL datasets (greater overall sample size, even if performed only for blood cells)
Other functional molecular QTLs such as methylation QTLs, histone QTLs, DNA accessibility QTLs, etc.

The largest eQTL datasets I can think of is probably eQTLGen. This is a meta-analysis of many studies using >30,000 samples from whole blood. This might be less ideal if you want cell-type-specific QTL information.

For other molecular QTLs, I do not think there are large-scale datasets yet.


AttributeError: module 'scater' has no attribute 'normalize' I am using rpy2 to run normalization using scran package. After calling computeSumFactors I want to do the actual normalization calling normalize() function from scater package:

def normalize(counts):
utils = importr('utils')
utils.chooseCRANmirror(ind=1)

packnames = ('Matrix', 'limSolve')
names_to_install = [x for x in packnames if not isinstalled(x)]
if len(names_to_install) > 0:
    utils.install_packages(StrVector(names_to_install))

base = importr('base')
base.source('http://www.bioconductor.org/biocLite.R')
biocInstaller = importr('BiocInstaller')

biocInstaller.biocLite('SingleCellExperiment', suppressUpdates=True)
biocInstaller.biocLite('scran', suppressUpdates=True)
biocInstaller.biocLite('scater', suppressUpdates=True)

sceLib = importr("SingleCellExperiment")
scranLib = importr("scran")
scaterLib = importr("scater")
named_list = robjects.r.list(counts=counts)
sce = sceLib.SingleCellExperiment(named_list)

try:
    clusters = scranLib.quickCluster(sce)
    sce = scaterLib.computeSumFactors(sce, cluster=clusters, positive=True)
except rpy2.rinterface.RRuntimeError:
    print("fewer cells than the minimum cluster size, so no quickCluster is run")
    sce = scranLib.computeSumFactors(sce, positive=True)

----> sce = scaterLib.normalize(sce, logExprsOffset=1)

Everything runs fine except the last line which gives the error:

  AttributeError: module 'scater' has no attribute 'normalize'

This is strange since that is what I am supposed to do actually according to docs.

I suspect that maybe I am loading a wrong scran version? I do not know though neither how to check the version number, nor how to install a specific one. Any suggestions would be greatly appreciated.

I tried the following:

sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)

And it is working. I am afraid though that it is a wrong way to go.

I tried to get the normalized matrix but it retrieved the initial raw counts one:

sce = scaterLib.normalizeSCE(sce, log_exprs_offset=1)
mtx_norm_r = robjects.r.assay(sce)
mtx_norm = np.array(mtx_norm_r)

So, something is definitely wrong here. Either I am retrieving the normalized matrix in the wrong way or the normalization was not run correctly. I tried the following method for retrieval but no luck:

mtx_norm_r = sceLib.normcounts(sce)

  RRuntimeError: Error in assay(object, i = exprs_values) :
    'assay(, i="character", ...)' invalid subscript 'i'
  'normcounts' not in names(assays())

I found a way to get the version numbers using getNamespaceVersion() function:

robjects.r.getNamespaceVersion('scater')

The versions are:

scater - 1.8.4
scran - 1.8.4
SingleCellExperiment - 1.2.0

They are the newest as I understand.

I was able to get to the normalized values using logcounts function in the following way:

mtx_norm_r = sceLib.logcounts(sce)
mtx_norm = np.array(mtx_norm_r)

Since I used log_exprs_offset=1 the field that I needed to use is logcounts, otherwise it should be normcounts which I found here.


How do you currently track bugs in python using freely available software/websites? I am starting a project using Python to track, archive, assign and manage bioinformatics bugs. This software will be able to be integrated into project management software or stand alone.

How do you currently track bugs/assign bioinformatics bugs using freely available bioinformatics software/websites?
What are features that you like/dislike? What are features that you would like to see in a bioinformatics bug tracker?

I am a Biology student with basic/core Python skills.

Use IPython and activate the trace in the Python debugger as follows
%pdb 1

or
%pdb on

to turn this off %pdb 0 or %pdf off
I do recommend using a formal IDE such as PyCharm or Visual Studio Code. The professional version of PyCharm is highly rated and within that use breakpoints . In Python3.7 there is supposed to be breakpoint() so you may not need an IDE, i.e. you set the breakpoint in your code using breakpoint(). In older versions of Python you can import this function via, import pdb; pdb.set_trace()
Finally check out PDB++ which is found here
Its important to note that experienced programmers often debug by simply inserting print() throughout their code and don't bother with an IDE


A program "PSICOV" "stopped" from a few days on Ubuntu 18.04? I am running dmpfold, "seq2maps.sh" and from a few days the program is at the "PSICOV" stage:

psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

"ps aux" gives me this:

daniel     26693  490  2.5 2546588 2514320 ?     Rl   سبت20 25116:58 psicov -z 15 -o -d 0.03 tail-fiber-sequence.aln

Any idea on what could be the matter? Yes, I am runinng "rosetta" too on my 32 CPUs, so all my CPUs are covered to max 100%; howver, RAM is 93Gb and only 32 GB has been used. Thanks

In summary just be patient. psicov doesn't  parrallelise so a tenfold increase in CPU will not hurry it along.

Also psicov is not RAM intensive given the machine you are describing, so you haven't hit a RAM bottleneck.

I don't know the ps output format you are using I'm use to seeing %CPU and %MEM from ps -aux on Susie and OSX. You should also try,

top -c -u username

or for OSX top -u

You should see the psicov happily chuntering away consuming ~100% of a single process. There is a UNIX/Linux guru specialist on BioSE who could almost certainly provide better diagnostics.


3D PCA group labelling I would like to make a 3D PCA but not sure how to label group wise which i can do for 2D PCA

library(ggplot2)
test<- read.csv("NEW_RBP/TF_DISEASE.txt",header = T,row.names = 1,sep = '\t')

Sample = c(rep("HSC",4),rep("Blast",11),rep("LSC",8))
Stage=c(rep("HSC",4),rep("Blast",11),rep("LSC",8))

df_final=data.frame(Sample,Stage)
df_final

pca_data=prcomp(t(test), center=TRUE, scale=TRUE)
pca_data_perc=round(100*pca_data$sdev^2/sum(pca_data$sdev^2),1)

df_pca_data=data.frame(PC1 = pca_data$x[,1],
                       PC2 = pca_data$x[,2], sample = colnames(test))
myColors <- c("red", "black", "blue2")
sample <- as.vector(df_pca_data$sample)
ggplot(df_pca_data, aes(PC1,PC2, colour = factor(Stage)))+#,label=rownames(t(test)))+
  geom_point(size=30)+
  geom_text(aes(label=sample),color="white",size=6,angle =0,parse = TRUE,  family="Bookman", fontface="bold")+
  labs(x=paste0("PC1 (",pca_data_perc[1],")"), y=paste0("PC2 (",pca_data_perc[2],")"))+
  theme_minimal(base_size=20) +
  theme(axis.text.x=element_text(size=rel(3), angle=90))+
  theme(axis.text.y=element_text(size=rel(3), angle=90))+
  theme(axis.title.x = element_text(colour="grey20",size=55,angle=0,hjust=.5,vjust=0,face="plain"))+
  theme(axis.title.y = element_text(colour="grey20",size=55,angle=90,hjust=.5,vjust=.5,face="plain"))+
  scale_color_manual(values=myColors)

The above code give me the figure as i have attached

For 3D PCA I tried this

scores = as.data.frame(pca_data$x)
head(scores)

plot3d(scores[,1:3], col=c(1:4), size=20, type='p',
       xlim = c(-50,50), ylim=c(-50,50), zlim=c(-50,50))
text3d(scores[,1]+2, scores[,2]+1, scores[,3]+1,
       texts=c(rownames(scores)), cex= 0.7, pos=3)

I get something like this as in case of 2D PCA i can label lets say all of the HSC sample as black and so on..Im not sure how to do that for 3D as I tried using factor i couldn;t do

Any help or suggestion would be highly appreciated


Are there genomics databases for neurological and psychiatric disorders? Are there databases that relate genetic disorders with neurological disorders?

I am thinking of something like COSMICS for cancer? AFAIK cancer usually develops due to a combination of 5-10 mutations.

Also if you can comment how genomics/bioinformatics in neurological contexts is compared to the how bioinformaticians study cancer?

Related: https://www.nimhgenetics.org/resources/psychencode

Try with PsyGeNet:

  PsyGeNET (Psychiatric disorders Gene association NETwork) is a resource for the exploratory analysis of psychiatric diseases and their associated genes

It is based on literature mining and curated by experts.

Here is the publication if you want more details:

PsyGeNET: a knowledge platform on psychiatric disorders and their genes

and they even have a Bioconductor package for programatic queries:

psygenet2r: a R/Bioconductor package for the analysis of psychiatric disease genes


Generate SNP/indel annotation in Arabidopsis I have a genome-wide list of germline SNPs and short indels for Arabidopsis thaliana, which I generated with Varscan. Regardless of the tool used to generate them, I would like to annotate them, i.e. knowing which ones can cause an aminoacid change or an early stop, using the default Arabidopsis thaliana Columbia 0 cultivar, for which I have both the sequence (FASTA from TAIR10) and the updated annotation (GFF from Araport).

Chrom   Position    Ref Var
Chr1    626503  G   T
Chr1    926694  C   T
Chr1    5280350 C   A
Chr1    5699993 C   A
Chr1    7004559 G   A
Chr1    8325810 C   T
Chr1    9371723 T   G

What I want to do is similar to what Annovar does, but unfortunately Annovar does not support Arabidopsis. I was thinking of an already existing R pipeline that takes in a genome, an annotation, a SNP/indel list and boom, annotation. But I couldn't find any, except maybe snpEffect. Any tips?

The Ensembl Variant Effect Predictor will annotate your A. thaliana variants:
http://plants.ensembl.org/Oryza_sativa/Tools/VEP

The web based tool will use the reference genome assembly and gene annotation used in the Ensembl browser:
http://plants.ensembl.org/Arabidopsis_thaliana/Info/Annotation/#assembly

However, if you use the command-line version of VEP, you can provide your own reference genome assembly and gene annotation files to be used as a reference.

You can find all of the documentation here:
http://www.ensembl.org/info/index.html

Best wishes

Ben
Ensembl Helpdesk


Filter bam using SNP list in bed format with minimum mapping quality and base quality I have a bam file and a bed file that defines a list of SNPs. I would like to filter the bam file to contain only those reads with a minimum mapping quality that overlap at least one SNP with a minimum base quality.

Samtools seems to almost solve this problem, but not quite:

samtools depth -b snp_list -q 20 -Q 30 example.bam allows me to count the number of alignments meeting my desired criteria, but not to produce a bam with these alignments.

samtools view -L snp_bed -q 30 example.bam allows me to filter the alignments overlapping any SNP with minimum mapping quality, but not with a minimum base quality at the SNP sites. There is no -Q option for samtools view, and the -q option is for mapping quality

Can this be achieved in samtools? Or is there another preexisting tool that can?

Thanks Devon ;-)
using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html

a bed file 'input.bed'

and the following snippet:

private IntervalTreeMap<Boolean> treeMap=null;
private final int BASE_QUALITY= 20;
@Override
public Object apply(final SAMRecord record) {
if(record.getReadUnmappedFlag() || record.getCigar()==null) return null;
if(treeMap==null) {
    treeMap = new  IntervalTreeMap<Boolean>();
    final String bedStr;
    try {
        bedStr = IOUtil.slurp(new java.io.File("input.bed"));
        }
    catch(java.io.IOException err)
        {
        throw new RuntimeIOException(err);
        }
    Arrays.stream(bedStr.split("[\n]")).
        filter(L->!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L->L.split("\t")).
        map(A->new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I->treeMap.put(I,true));
    }

byte quals[]= record.getBaseQualities();
int refpos= record.getStart();
for(CigarElement ce:record.getCigar()) {
    CigarOperator op= ce.getOperator();
    if(op.consumesReferenceBases())
        {
        if(op.consumesReadBases()) {
            for(int i=0;i< ce.getLength();i++)
                {
                Interval t = new Interval(record.getContig(),refpos+i,refpos+i);
                if(!treeMap.containsOverlapping(t)) continue;
                int readpos = record.getReadPositionAtReferencePosition(refpos+i);
                if(readpos==0 ) continue;
                readpos--;
                if(readpos>=quals.length)  continue;
                if(quals[readpos]>BASE_QUALITY) return true;
                }
            }

        refpos+=ce.getLength();
        }

    }
return null;
}

usage:

samtools view -bu -L input.bed input.bam |\
   java -jar dist/samjdk.jar --body -f snippet.code


Installing DESeq2 in Ubuntu I am trying to install DESeq2 in my Ubuntu with R version 3.5.1. According to the package repository in Bioconductor the version should be 3.5.

> R.version

platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

then I attempt to install DESeq2, following the instructions in Bioconductor:

> if (!requireNamespace("BiocManager", quietly = TRUE))
       install.packages("BiocManager")
   BiocManager::install("DESeq2", version = "3.8"

and at the end I get:

Warning messages:
1: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘RCurl’ had non-zero exit status
2: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomeInfoDb’ had non-zero exit status
3: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘annotate’ had non-zero exit status
4: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomicRanges’ had non-zero exit status
5: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘genefilter’ had non-zero exit status
6: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘geneplotter’ had non-zero exit status
7: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘SummarizedExperiment’ had non-zero exit status
8: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘DESeq2’ had non-zero exit status

I have visited $2*10^{9999}$ discussions on this issue and none of them helped me solving the problem. In addition, I didn't have any problem installing other packages like limma and edgeR, so I assume that this might be specific to DESeq2 or one of it's dependencies.

Thank you for the help!

First I installed this libraries:

sudo apt-get install libcurl4-openssl-dev libxml2-dev

I then installed RCurl:

BiocManager::install("RCurl")

and lastly:

BiocManager::install("DESeq2", version = "3.8")


Human Herpes Virus, fusion motif in protein vs. SARS-CoV-2? Is it possible that the fusion with cell propriety of Human Herpes virus, I see X14112.1 in GenBank is due to protein UL36 (tegument of virus) starting in linear position 2685 ???
H R R A R S L A

there are similarities with SarsCov2 virus protein S in name MN908947.3 in GenBank position 681.
P R R A R S V A

It is strange but the now (21-12-2020) the new Sars2 Great Britain strain is more near to that UL36 piece, it seems already:
H R R A R S V A

Herpes virus protein UL36 and SarsCov2 virus protein S[1200 position], have  too
L Q E L G K

i find in internet they call it "Glutaminase kidney isoform, mitochondrial"

I can answer that ... I'm afraid no. What you are asking is whether there is any similarity in their cellular receptor binding or "fusion loops" because the motifs look a bit similar.
These are very different viruses and you need to appreciate just how different they are,
Human herpes is a double stranded DNA virus resulting in a chronic infection. It has a vast genome of >150kb, which includes cytokine mimics and proteins that we have no idea what they do. The protein you are referring to is likely a surface antigen in a virus covered in different surface antigens. However, this protein is involved in blocking the beta inferon response is considered a protease - not fusion loop and appears to have some sort of complex replication function involving the capsid. Finally, there is also the question of which human herpes 1 to 8. Here you are referring to HHV1 but for example HHV3 is chickenpox (okay its an acute not chronic infection) and HHV5 is cytomegliovirus (congenitcal infection).
COVID-19 is a simple virus of a 30Kb genome of a single stranded positive sense RNA, causing an acute infection. I agree the Spike protein is much bigger than most RNA viruses, but the receptor-binding is tightly defined with ACE-2. It forms a classic RNA icosahedral virus, singly covered in Spike protein. I have no doubt it will go through conformational change like other RNA viruses in receptor binding but the protein involved is singularly involved in cell entry - thats its job.
Summary
There are large differences between HHV1 and COVID-19 for,

Viruses; DNA double stranded vs RNA single stranded
Infections; chronic vs. acute
Transmission; vertical/integration? direct contact vs. respiratory
Ultimately very different proteins
Finally very different functions

At best this would be heavily  scrutinised as a virological assumption.
Bioformatics If you really wanted to follow this up you could compare the structures around these motifs, because COVID-19 is close to SARS and assess this possibility ... but you would need to be careful about the conformational changes Spike goes through in receptor-binding.
Advice before investing large amounts of time in a bioinformatics question, it is well-worth understanding the biological rationale for the question. In this instance that rationale appears to be lacking. If you found a comparison for example with gamma-coronaviruses, which was absent from many beta-coronaviruses that would be interesting. Here there is simply too much possibility of a chance similarity and given the size HHV1 that cannot be easily overlooked. However it is good that you are thinking and questioning, if you apply the same techniques, e.g. psi-blast, protein-blast into different questions ... that can work.


How to identify gene expression signatures from gene expression data? I have TCGA gene expression data. I'm interested in identifying gene expression signatures using the data.

I would like to know whether there are any tools or R packages for identifying gene signatures.

How are gene signatures different from GSEA?

There is no golden/standard way to define a gene signature but they are completely different from a gene set enrichment analysis (GSEA). I will start with how to obtain a gene signature:

Usually this kind of signatures are defined by comparing a group against another ie doing a differential gene expression analysis. One then selects the most up or down regulated genes and defines them as the signature of the group of cells. That's how the gene sets of the MSigDB which are not from (or derived from) pathway databases or gene ontologies are defined. There is no gold threshold of a fold change or a specific number of genes, it (highly) depends on the experiment you run and how good are the isolation of the cells.
Tools: DESeq2, limma, edgeR,...

As suggested by gringer, another way is to find which genes are expressed together (independently of how do they relate to other groups). This method relays on having a good enough number of samples from the same condition to be able to identify which genes are expressed the same way in those cells. If these genes are not expressed that way in other cell lines, then you have your signature for this group. Again, there isn't a golden method, different types of co-expression networks can be build and depends on how specific do you want to be or which is the data provided.
Tools: WGCNA, diffcoexp, TReNA,...

This is different from a gene set enrichment analysis1, which consists on evaluating if a group of genes (usually a gene set, a pathway or genes from a specific gene ontology) are not randomly sorted in a sorted list of genes (usually this sorted list is the fold-change of a differential analysis).

There are several methods developed depending on the null hypothesis, the type of sub-method and what do they take into account...
Tools: GSEA, PLAGUE, GSVA, limma, fgsea,...

TLDR: To perform a gene set enrichment you don't need a gene signature but a gene set (that can be a gene signature or a pathway or a random group).

1: Note that gene set enrichment analyisis is different from gene enrichment analysis. For a good review on the methods I recommend this article


Can blat use more than one core/CPU to speed up the alignment? I am using BLAT to align two versions of the genome of C. elegans. I can see in the Activity Monitor of my Mac Book Pro High Sierra that blat is using 100% of a CPU. However, is this programme able to use more than one core / CPU to speed up the alignment?

I cannot find such option in the help page of the tool here.

I installed blat using conda install -c bioconda blat.

I am using the latest version of blat, i.e. 36.

BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.

EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based pipeline. So far there are no converters for minimap2's SAM/PAF or MUMmer's delta formats. With the UCSC pipeline, you have to write a converter by yourself. An alternative approach is to use minimap2's liftover. The minimap2 cookbook gives an example:

minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -

Note that this pipeline is not as mature as UCSC liftOver and is not as scalable.


Getting species from ensembl sequence id I have downloaded some orthologues from ensembl website in fasta format. The problem is that it uses the ensembl sequence id as name for the sequences. I want to know the species name or any other information about the sequence (chromosome, location...) For example for this sequence id: ENSMODP00000018933. How can I do that? I have been looking at the REST API but the most I can find is an API call for getting the sequence passing the id as parameter.
https://rest.ensembl.org/documentation/info/sequence_id
I have also tried this API call but does not work for sequence ids:
https://rest.ensembl.org/documentation/info/xref_id

The lookup/id endpoint will get it for you. Of you can just look up the three letter species code, in this case MOD, on the list in the documentation.