Create new range based from existing range (python) I have a pyranges object that I am looking to iterate over. For each range, I would like to "resize" it to its centre (so start+end/2) and save that into a separate pyranges object. There will be additional processing associated based on what ranges already exist nearby, so I can't just use a resize function, I have to iterate over the object.
First, how would I go about manually creating a new pyranges object, then inserting one by one new ranges into it?
Second, when I try to iterate over my current pyranges object (using this line for myspot in pr.itergrs(hotspots, strand=False) :), it gives me this error: AttributeError: 'tuple' object has no attribute 'stranded'. This error occurs whether or not I have strands in my pyranges, and whether or not I specify strand=False. Why does this occur?

The data I am working with: a pyranges covering hotspots in the human genome (all chr1 for the moment)
+--------------+-----------+-----------+-----------+
| Chromosome   | Start     | End       | Score     |
| (category)   | (int32)   | (int32)   | (int64)   |
|--------------+-----------+-----------+-----------|
| chr1         | 154421820 | 154421821 | 26413     |
| chr1         | 30749469  | 30749470  | 14833     |
| chr1         | 55259503  | 55259504  | 9900      |
| chr1         | 14985097  | 14985098  | 7735      |
| ...          | ...       | ...       | ...       |
| chr1         | 247975751 | 247976497 | 1         |
| chr1         | 247977599 | 247981391 | 1         |
| chr1         | 248018012 | 248019454 | 1         |
| chr1         | 248177298 | 248178490 | 1         |
+--------------+-----------+-----------+-----------+

The hotspots are sorted by their score, and I'd like to go through and for each one, resize it to its centre, keeping the score, but if there is already a hotspot within a certain range (say +/-5kb), to discard it instead, so that I end up with clear chunks between hotspots. This is the "pseudocode" I was planning to do to implement this:
add centre of hotspot 1 to a new pyranges object
create ranges with score = -1 for 5kb flanking each side
for myhotspot in hotspots:
  resize to centre
  if not overlap with new pyranges object
    add centre
    create flanking regions
discard ranges with score<0

This may not be the most elegant solution, but I'll never have more than ~2000 ranges to process so it should work fine. Is this implementable with pyranges?

Answer to main question:
import pyranges as pr
assert int(pr.__version__.split(".")[2]) >= 93, "pip install pyranges==0.0.93"
import numpy as np
np.random.seed(42 * 10)

# create large df to test on
gr = pr.random(int(1e5), length=1, chromsizes={"chr1": 249250621})
gr.Score = np.random.randint(250, size=len(gr))

result = gr.max_disjoint(slack=5000)

# small example that is easy to validate
gr2 = pr.from_string("""Chromosome Start End Score
chr1 1 2 5
chr1 4999 5000 10
chr1 10001 10002 1
chr1 15000 15001 2
chr1 20100 20101 5
""")

gr2.max_disjoint(slack=5000)

# +--------------+-----------+-----------+-----------+
# | Chromosome   |     Start |       End |     Score |
# | (category)   |   (int32) |   (int32) |   (int64) |
# |--------------+-----------+-----------+-----------|
# | chr1         |         1 |         2 |         5 |
# | chr1         |     10001 |     10002 |         1 |
# | chr1         |     20100 |     20101 |         5 |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 3 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Could you please post an example of what data you have and what you want to end up with?

Inserting row after row is going to be slow, whether it is pandas or pyranges. It is best to collect all the data and then create a new dataframe or PyRanges.

pr.itergrs is for iterating over multiple pyranges at the same time. It takes a list of pyranges while you seem to have a collection of tuples. The function iterates over every chromosome/strand pair in the list and for it extracts the data from every pyranges. See the code example at the end.

To iterate over a pyranges you can do for k, df in gr: # ....
Anyways, to do what you want:
import pyranges as pr
import numpy as np

f1 = pr.data.f1()
f1.Start = ((f1.Start + f1.End)/2).astype(np.int32) # int is 64 bit
f1.End = f1.Start + 1

I understand I have not done the addititonal processing you speak of, but if you can post an example I'm sure I'd be able to help you :)

Code example pr.itergrs:
import pyranges as pr
f1 = pr.data.f1()
f2 = pr.data.f2()
for k, dfs in pr.itergrs([f1, f2], keys=True):
    print(k)
     for df in dfs:
         print("df:")
         print(df)

('chr1', '+')
df:
  Chromosome  Start  End       Name  Score Strand
0       chr1      3    6  interval1      0      +
2       chr1      8    9  interval3      0      +
df:
  Chromosome  Start  End Name  Score Strand
0       chr1      1    2    a      0      +
('chr1', '-')
df:
  Chromosome  Start  End       Name  Score Strand
1       chr1      5    7  interval2      0      -
df:
  Chromosome  Start  End Name  Score Strand
1       chr1      6    7    b      0      -


Is BWT based aligner suitable for any types of alignment task? Burrows wheeler transformation based aligner like BWA or bowtie seems a standard alignment tool used many area. I was just wondering if there is a kind of alignment task in which BWT algorithm is not suitable.

BWT based algorithm is just a component of a full-pledge aligner. To that end, it is applicable to most alignment tasks. However, for long-read alignment and metagenomes, we don't often use BWT based algorithms as other methods are faster.


How to dump genes from GenBank in GFF3 format? This question has also been asked on BioStars

If I look at this record in GenBank I see about 6k genes:
https://www.ncbi.nlm.nih.gov/nuccore/CM000760?report=gbwithparts

I'd really like to be able to dump those genes in GFF3 format, but I'm guessing I can't do that in a single step, so...

How to dump all genes from a GenBank record in 'whatever' format?
How to convert 'whatever' format to GFF3?

Turns out you can just grab the GFF3 from the NCBI's FTP site!

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/003/195/GCA_000003195.3_Sorghum_bicolor_NCBIv3/

See:

https://www.ncbi.nlm.nih.gov/genome/doc/ftpfaq/#files

Or "Access the data" on the right here:

https://www.ncbi.nlm.nih.gov/assembly/GCF_000003195.3


Existing tool for converting gff3 to genbank (gbk) I want to convert my gff3 annotation files to genbank format for use in Mauve. I found the seqret tool here https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/ which can perform this task, but my files (bacterial genomes) are too big.

Can anyone suggest a simple method to make this conversion?

I believe seqret is already the simplest approach. If you install the emboss command line utility, genome (file) sizes will not be a problem. You can even install it with conda, see https://anaconda.org/bioconda/emboss.

After that, you can use it as the example:

seqret -sequence {genome file} -feature -fformat gff -fopenfile {gff file} -osformat genbank -osname_outseq {output prefix} -ofdirectory_outseq gbk_file -auto

Hope it helps


Install DESeq2 through anaconda I tried running conda install -c bioconda bioconductor-deseq2 in a conda environment, but when I run R console and try to import DESeq2 I am getting the error message:

  Error: package or namespace load failed for ‘GenomeInfoDb’ in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
   there is no package called ‘GenomeInfoDbData’
  Error: package ‘GenomeInfoDb’ could not be loaded

When I try installing GenomeInfoDbData and GenomeInfoDb using conda, it says that the packages were already installed. If I run install.packages(GenomeInfoDb) from R console, it is telling me that the package is not available for R version 3.5.1.

I want to import DESeq2 in Rstudio though and it is giving me a different error when I try importing DESeq2:

  Error in library("DESeq2") : there is no package called ‘DESeq2’

I do not get why Rstudio does not see that I alrerady installed DESeq2 with conda. If I try running biocLite("DESeq2") from Rstudio I get the following error:

  Error in dyn.load(file, DLLpath = DLLpath, ...) :
    unable to load shared object '/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so':
    dlopen(/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so, 6): Symbol not found: _libintl_dgettext
    Referenced from: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so
    Expected in: /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/lib/libR.dylib
   in /Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/cluster/libs/cluster.so
  ERROR: lazy loading failed for package ‘Hmisc’
  * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/Hmisc’
  ERROR: dependency ‘Hmisc’ is not available for package ‘DESeq2’
  * removing ‘/Users/nikitavlasenko/miniconda3/envs/renv/lib/R/library/DESeq2’

So, supposedly the issue is with Hmisc. I installed the package successfully with conda, but Rstudio is apparently does not know about it. Running install.packages("Hmisc") gives the first part of the error with missing cluster.so library.

Any suggestions would be greatly appreciated.

Ultimately my colleague helped me to solve the issue by following the steps:

Created environment: conda create --name myenv
Activated it: source activate myenv
Edited .condarc to set the priority of the channels to be the following (top to bottom): conda-forge, bioconda, r, default
Checked that the channels are set in the correct order: conda config --get channels
conda install r=3.5.1
conda install r-essentials
conda install rstudio
Installation of all of the bioconductor packages through bioconda and never install anything from inside R or Rstudio (I guess that is how it was broken)

Then launching rstudio from within the environment. Maybe it is not a kind of quick fix and it could have been solved easier, but I was not able to.


Convert paired-end BAM into a single-end BAM and keep all the reads I want to call peaks using MACS2 on my paired-end ATAC-seq BAM and I am interested in cutting sites. So I make use of the -BAM option (opposed to the -BAMPE) option. However this throws away all the mates, and all the reads of which the mate did not align (that's more than 50%!!!). My idea was to align paired-end, and then remove all the information from the BAM that it is paired-end. So my question is how do I:

Convert paired-end BAM into a single-end BAM and keep all the reads

I made a start based on the SAM-format specification , but when I run the result of this script through markduplicates I get the error:

ERROR: Record 1, Read name SRR6245286.60336908., MRNM should not be set for unpaired read.

I can not find what this MRNM is...

This is my code I used to generate the faulty 'single-end' bam, but I am open to suggestions (e.g. bash, preferably no R).

import pysam

paired =         1
proper_pair =    2
mate_unmapped =  8
mate_reverse =   32
first_in_pair =  64
second_in_pair = 128

bam_in = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7.samtools-coordinate.bam", "rb")
bam_out = pysam.AlignmentFile("/vol/atacpipeline/bwa/GSM2837489-danRer7_SE.samtools-coordinate.bam", "wb", template=bam_in)
for line in bam_in:
    if line.flag & second_in_pair:
        line.pos += line.template_length
    line.next_reference_id = 0
    line.next_reference_start = 0
    line.flag &= ~(paired + proper_pair + mate_unmapped + mate_reverse + first_in_pair + second_in_pair)

    bam_out.write(line)

Like Devon pointed out, most likely you should sort out whether the files have been marked for duplicates correctly.

You can also use samtools rmdup too

samtools rmdup <bamfile> <output>

Also, I find it very odd that 50% of the reads are lost when you keep only properly paired reads.

I would suggest removing the duplicates and then converting the bam file to a bed file. Macs2 can take in the bed format. Any extremely crude python script that converts a bam into a bed:

import pysam
samfile = pysam.AlignmentFile("test.bam", "rb")
open("test.bed", 'w').close()
f = open("test.bed", "a")
STRAND = ["+","-"]
for read in samfile.fetch():
    STR = STRAND[int(read.is_reverse)]
    BED =[read.reference_name,read.pos,read.reference_end,".",read.mapq,STR]
    f.write('\t'.join([str(i) for i in BED]))
    f.write('\n')
f.close()

You can also try the conversion with bedtools or definitely refine that script above with pandas etc.


Source for whole genome comparisons at NCBI Genomes The NCBI Genomes database has these dendrograms for (presumably) whole genome comparisons for certain species, e.g. Pseudonomas aeruginosa  or Escherichia coli.

How were these comparisons done? Someone knows the source / paper?

I inquired into the details of the dendrograms after becoming frustrated with the lack of information. As with Ensembl, I'm sure that the folks at NCBI have a standardized pipeline that they run the sequences through to generate these dendrograms as no specific source deserving attribution was mentioned:

The tree is based on a pairwise, BLAST comparison of chromosome sequences (an assembly is used for genomes with 2 or more chromosomes). BLAST identity is normalized on sequence length. The distance map is generated by a Neighbor Joining algorithm, where the distance is the blast score.
Dr. Wayne Matten, NCBI


How many reads do I need for hybrid assembly I have Illumina and PacBio reads and I would like to use dbg2olc for hybrid assembly.

Part of dbg2olc is SelectLongestReads which select reads that sum to total_length of bases.

How many bases do I need for a 1.8G allotetraploid to remain?

I have never used dbg2olc. However other OLC assemblers managed to do a solid assembly from 20x of corrected long reads. But obviously, the more the better.

In you case of allotetraploid, you need to figure out if the heterozygosity is sufficient to separate haplotypes of not and adjust parameters of the assembler. For instance in bdelloid rotifers the tetraploid genome has an AA'BB' structure, all the assemblies were able to separate A and B, but not A and A' and B and B'. I guess an allotetraploid could have a similar structure.

It would be good to know what you are dealing with and what is your goal in the assembly. Do you want a single haplotype?

-- edit --

You can check the genome structure using short reads by smudgeplot, a tool I made (yes, it's a shameless plug).


Relationship between the instantaneous rate of change of a continuous time Markov model and the rate of a global clock model? If continuous-time Markov models contain a parameter q that denotes the instantaneous rate of change and a global/strict clock model is the overall rate of evolution, how are these two parameters related? That is, for a given substitution model (e.g., JC69, F81, GTR), is there a deterministic relationship between the rate parameter(s) of these substitution models and the rate of the clock model? For instance, should the individual instantaneous rates of change (e.g., A -> C, T -> C, etc.) sum up to the global clock rate? If this is not the case, how should one understand the relationship between the rate of a continuous-time Markov model and the overall evolutionary rate?

Essentially you are comparing a "local clock" with a global clock and whether the instance rate relates to the global it depends on the model and in MCMC clocks in most cases the answer is no. A global clock is a very specific and very limiting model (below). Everything below is a comparison between taxa not between alignment positions.

"Local clock" in the context of MCMC clocks is in inverted commas because they are using population models. I will use the nomenclature here for convenience to describe rate heterogeneity between branches, not within a branch (last example).

A "local clock" has mutation rate heterogeniety across the tree whilst a global clock is completely homogenious. You can't say how is a "local clock" related to global clock because there are numerous different models.

Global clocks never hold for biological data because there are always some taxa evolving at a faster rate than others. This is particularly true in infectious disease.

One historic model used in traditional local clocks (no inverted commas) was autocorrelation, where a local clock would relate to a global clock on a "clade by clade" basis. That is clades are correlated by their mutation rate, so some clades evolve faster and other clades slower. Statastically an autocorrelated clock always outperforms a global clock under a likelihood ratio test.

I remember autocorrelation clearly because I implemented a variant of this approach, so I needed to understand how they did it, rather than click a few buttons on a GUI, which is what happens in Beast.

Beast has at least 3 population models and off-hand I cannot remember the population scenarios. I do remember the population scenarios made sense and secondly they can give 10-fold differences in rate in the data set I looked at, You need to read on the background theory for each model and assess what is more biologically appropriate for a given context. One model (from dim memories) reminded me of a classic population growth equation, i.e. an exponential. Population growth and collapse in phylodynamics would be informative of the rate, but again you need to go to the underlying theory.

The molecular clock theory that the rate of evolution towards the tips of trees is much faster than at the deeper branches is the sort of model that could be implemented in a MCMC clock. In the latter case the relationship between a global clock and the rate heterogeneity clock (which is a better description in that model than "local clock") is quite different from global vs. "local clock".

The fundamental issue with clocks is exemplified in the 10-fold difference issue. The lesson that is important to me is very easy to get a "wrong answer". The lesson important to you is if two MCMC models can't agree the same rate how can both obtain the same answer as a global clock? I am personally a fan of maximum heterogeneity, but in this day and age they will ask you to perform a Beast model to support it.


R packages for data analyses of pooled CRISPR screens We are designing a CRISPR/Cas9 experiment and thinking of the down-stream data analyses.

Are there any R packages for analysing raw NGS read count data from pooled genetic screens using CRIPSR/Cas9 to disrupt gene expression in a population of cells? I guess we will need to start with basics, i.e. sequences processing, data exploration, visualisation etc.

Assuming you mean CRISPR screens targeting many loci (for example using the GeCKO library) there is an R package here.

It uses a linear mixed effect model to compare guide counts before and after a selection step, allowing for multiple guides/gene and multiple replicates.
It can also do the initial read alignment AFAIK and the author should be quite responsive in case you have questions.

Another option could be the python program MAGeCK


How to select specific lines and specific columns Please I'm a beginner in using R and I work on database.
I have  a database that contain lot of columns and lines, like this for example
Column1   Column2
5            T1
1            T2
1            T3
1            T1
2            T1
1            T1
2            T1
1            T1

I want to extract just lines that have the value 1 in column 1 and Value T1 in Column 2
So the result I want is
Column1   Column2
1            T1
1            T1
1            T1

Can you tell me please how I can do that
Thank you very much

You can simply use subset function in R and use the code given below,
    new_dataframe <-
    subset(name_of_your_dataframe,Column1 == 1 &&
    Column2 ==  "T1")

Please keep in mind that R is case sensitive and you should be careful while using the variable names.


Error downloading GSE dataset I've installed GEOquery with
>install.packages("BiocManager")
>BiocManager::install("GEOquery")
>library(GEOquery).

then I ran the command
>GSE68808 <- getGEO("GSE68808")

but it resulted in..
>Found 1 file(s)
GSE68808_series_matrix.txt.gz
trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE68nnn/GSE68808/matrix/GSE68808_series_matrix.txt.gz'
Content type 'application/x-gzip' length 4138819 bytes (3.9 MB)
downloaded 3.9 MB

error: 'C:\Users\(myname)\AppData\Local\Temp\RtmpUp5Sv7\file8e844a7a59bc' does not exist.

ㅜㅜ
How can I solve this? thanks.

There seems to be some problem with your "platform specific tempdir" that is used by getGEO's destdir parameter. See ?getGEO for more info on the destdir parameter.
Try
GSE68808 <- getGEO("GSE68808", destdir = getwd())


Alternative splicing positions on protein I have mass-spectrometry data (Human cellsanalyzed and annotated by MaxQuant).
I want to find which peptides are covering positions of former splice junctions (on the mRNA sequence), based on protein id and peptide positions.

Where could I find an appropriate data? I want to know, based on peptide-sequence, annotaion and position on protein, whether this peptide is a single exon, or if it contained areas of former splice junctions.

Itai

If you are considering the human genome as "data" (not mRNA), I would run your single peptide via tblastn which will take your protein(s) translate them to nucleotides and run a blast against Genbank. The splice sites will be immediately apparent because they force indels in your query. If there are no indels there are no introns (splice sites). If you want to process several proteins, one solution is to use eblast, a command line version of tblastn  (Anaconda package). You can confirm this against a full genome annotation (note, I don't work on humans).

Alot of work has been done on splice sites, a friend worked on this, but I think you just want to focus on a restricted number of target peptides in your mass spec. output.

mRNA should be avoided IMO if you are looking for splice sites, i.e. intron boundaries.


RNA-seq gene/transcript read counts database for Mouse Is there any RNAseq gene/transcript read count database for the mouse? I already know about ARCHS4, looking for some other source.

Thank you

EBI has an expression atlas, including mouse: https://www.ebi.ac.uk/gxa/experiments?species=mus+musculus

If you click on one of the individual experiments you can download tables of TPM or RPKM per gene.

There are currently about 1000 individual experiments in different tissues etc.


Practical use of RNA structure I have spent some time studying RNA structures. While there is a range of interesting phenomena and functions that certainly deserve understanding from the scientific point of view, I have never encountered a practical use of this knowledge (in terms of disease diagnostics, medications, etc.) Have you? If a practical use of RNA structures is still in the future, what could that be?

Please look up flavivirus 'double loops' as you described them previously (post for "Coronavirus RNA') and associated RNA secondary structure anomalies for dengue virus and associated vaccine (Butantan) and the yellow fever virus and its vaccine (17D). If you are aware of "double loops", we must be aware of the association of RNA secondardy structure and attenuated strains, i.e. live vaccines. The point is that modification of the RNA secondary structure is a key strategy in flavivirus attenuation.


What's the most efficient file format for the storage of DNA sequences? I'd like to learn which format is most commonly used for storing the full human genome sequence (4 letters without a quality score) and why.

I assume that storing it in plain-text format would be very inefficient. I expect a binary format would be more appropriate (e.g. 2 bits per nucleotide).

Which format is most common in terms of space efficiency?

Genomes are commonly stored as either fasta files (.fa) or twoBit (.2bit) files. Fasta files store the entire sequence as text and are thus not particularly compressed.

twoBit files store each nucleotide in two bits and contain additional metadata that indicates where there's regions containing N (unknown) bases.

For more information, see the documentation on the twoBit format at the UCSC genome browser.

You can convert between twoBit and fasta format using the faToTwoBit and twoBitToFa utilities.

For the human genome, you can download it in either fasta or twoBit format here: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/bigZips/


Basecalling process carried out by full_1dsq_basecaller.py/.exe I am unable to understand what exactly the full_1dsq_basecaller.py/.exe script does. Is it giving the consensus sequence for linear and complementary strand, or it is just giving the reads that came from the complimentary strand?

If we are getting two fastq files as output which one we should use for alignment?

Thank you in advance.
I greatly appreciate any help

My understanding is that the 1D$^2$ process is carried out in three stages:

Carry out a 1D calling of reads (i.e. what the standard Albacore does at the moment)
Link reads that are similar to each other (in reverse-complement orientation), from the same channel, and have a very short pore transition time
Create a 2D consensus base call, using information from the two linked reads

The two outputs from the basecaller are the 1D calls from step 1, and the consensus calls from step 3. Ideally just the 2D consensus calls should be used for assembly (assuming it represents a random distribution of reads), but  if that gives too low coverage (or is not random enough), then the non-converted reads from the 1D calls could be included as well.

Because the linking is not a physical link, it's not a good idea to use 1D$^2$ for amplicon reads; the likelihood of two complementary reads from different templates appearing one after another is too high.


How to best detect the "peaks" in RNA-seq data that are not assigned to any gene? I encountered that many reads from single-cell RNA seq data were lost in the analysis because not assigned to any gene (genome: galgal6). I am trying to find an approach than could give me all the "peaks" in this data, that is to say regions where there is a high density of reads. Once I get there, I would extract "peaks" that are not assigned to any gene, and get the distance to the closest genes.

I am thinking of using first bedtools genomecov to get the "peaks", and then bedtools intersect to get all reads assigned to a gene. Finally, I extract reads that are not in the output of bedtools intersect to get a new file with all my reads from the peak, that are not assigned to any gene.

Does it make sense this way ? Can you see any other way of doing ?

Here are some tiny data:

> cat reads.bed
chr9    505479  505498
chr9    508014  508037
chr9    514603  514633
chr9    529519  529540
chr9    529519  529540
chr9    529519  529540

> cat tiny.galGal6.chrom.sizes
chr9    24153086

> bedtools genomecov -bg -split -i reads.bed -g tiny.galGal6.chrom.sizes
chr9    505479  505498  1
chr9    508014  508037  1
chr9    514603  514633  1
chr9    529519  529540  3

The last line would correspond to a peak. However, another difficulty is how should I define the thresholds ?

Any help would be more than welcome.

Basically what you're discovering is that there are unannotated expressed features, so your task isn't really finding peaks, but rather finding novel expressed transcripts. For that, you can use stringTie or similar programs. Ensure that the BAM file you give to stringTie have all of the cells together.


Feature annotation: RefSeq vs Ensembl vs Gencode, what's the difference? What are the actual differences between different annotation databases?

My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations. Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

To add to rightskewed answer:
While it is true that:
Gencode is an additive set of annotation (the manual one done by Havana and an automated one done by Ensembl),
the annotation (GTF) files are quite similar for a few exceptions involving the X chromosome and Y par and additional remarks in the Gencode file (see more at FAQ - Gencode).

What are the actual differences between different annotation databases?

They are a few differences, but the main one  for me (and it could be stupid) is
that Refseq is developed by the American NCBI and
the ENSEMBL is mainly developed by the European EMBL-EBI.
Often, labs or people will just start using what is the best known to them (because of a course or workshop) or because they start working with one of the databases with one specific tool and keep with it later.

My lab, for reasons still unknown to me, prefers Ensembl annotations (we're working with transcript/exon expression estimation), while some software ship with RefSeq annotations.

Your lab might be mostly European based people or they might also have read papers like the one from Frankish et al.  Comparison of GENCODE and RefSeq gene annotation and the impact of reference geneset on variant effect prediction. BMC Genomics 2015; 16(Suppl 8):S2 - DOI: 10.1186/1471-2164-16-S8-S2
From the Frankish et al. paper paper:

The GENCODE Comprehensive transcripts contain more exons, have greater genomic coverage and capture many more variants than RefSeq in both genome and exome datasets, while the GENCODE Basic set shows a higher degree of concordance with RefSeq and has fewer unique features.

As for:

Are there significant differences between them today, or are they, for all intents and purposes, interchangeable (e.g., are exon coordinates between RefSeq and Ensembl annotations interchangeable)?

No. I don't think they are great differences between them as that the global picture should stay the same (although you will see different results if you are interested in a small set of genes). However, they are not directly interchangeable. Particularly as there are many versions of Ensembl and Refseq based on different genome annotations (and those won't be interchangeable between themselves either in most cases).
However, you can easily translate most[1] of your Refseq IDs to ENSEMBL IDs and vice-versa with tools as http://www.ensembl.org/biomart/martview for example (there are devoted libraries/API as well like Biocondutor: biomaRt
[1] Most as sometimes, they might be annotated in one of the database but haven't (yet) an equivalent in the other.
EDIT
In fine, even if people tends to keep to what they are used to (and that the annotations are constantly expanded and corrected) depending on the research subject one might be interested in using one database over another:
From Zhao S, Zhang B. A comprehensive evaluation of ensembl, RefSeq, and UCSC annotations in the context of RNA-seq read mapping and gene quantification. BMC Genomics. 2015;16: 97. paper:

When choosing an annotation database, researchers should keep in mind that no database is perfect and some gene annotations might be inaccurate or entirely wrong. [..] Wu et al. [27] suggested that when conducting research that emphasizes reproducible and robust gene expression estimates, a less complex genome annotation, such as RefGene, might be preferred. When conducting more exploratory research, a more complex genome annotation, such as Ensembl, should be chosen.
[..]
[27] Wu P-Y, Phan JH, Wang MD. Assessing the impact of human genome annotation choice on RNA-seq expression estimates. BMC Bioinformatics. 2013;14(Suppl 11):S8. doi: 10.1186/1471-2105-14-S11-S8.


How does htslib/samtools access optional BAM fields? I am confused regarding how various software libraries deal with optional fields in a BAM:

Based upon the BAM specification, there are 11 mandatory fields to a BAM:

QNAME, FLAG, RNAME, POS, MAPQ, CIGAR, RNEXT, PNEXT, TLEN, SEQ, QUAL

The BAM files of today normally have many optional fields as well, which must be of the format TAG:TYPE:VALUE, whereby TAG is a two-character string that matches /[A-Za-z][A-Za-z0-9]/ and TYPE is a single case-sensitive letter which defines the format of VALUE. Each TAG can only appear once in one alignment line.

Certain BAMs may have no optional fields, while others have several.

My understanding is that optional tags may be in a different order from row to row, and then sometimes they simply do not exist.

My questions:

How does htslib/samtools (e.g. samtools view) access these optional fields, which differ in number from BAM to BAM (possible row to row)?

In htslib/sam.h, a BAM alignment struct is defined as follows:

/*! @typedef
 @abstract Structure for one alignment.
 @field  core       core information about the alignment
 @field  l_data     current length of bam1_t::data
 @field  m_data     maximum length of bam1_t::data
 @field  data       all variable-length data, concatenated; structure: qname-cigar-seq-qual-aux
typedef struct {
    bam1_core_t core;
    int l_data;
    uint32_t m_data;
    uint8_t *data;
#ifndef BAM_NO_ID
    uint64_t id;
#endif
} bam1_t;

and followed by several methods which access mandatory BAM fields, e.g.

/*! @function
 @abstract  Get the name of the query
 @param  b  pointer to an alignment
 @return    pointer to the name string, null terminated
 */
#define bam_get_qname(b) ((char*)(b)->data)

How does htslib access each of the optional fields, given this structure? Does it keep the order of the optional fields, or is this information unavailable?
When a user runs htslib or even samtools view, how does the BAM row "end", i.e. what data structure exists to show that the BAM row has no more optional fields, and that the algorithm should "move on" to the next BAM row.

Any clarification for how the source code works in this case is much appreciated.

How does htslib access each of the optional fields, given this structure?

htslib uses a macro, bam_get_aux, for this purpose. You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data. But you shouldn’t care about this, and just use the relevant exported functions for reading optional (“auxiliary”) data:

bam_aux_get returns a pointer to an optional record, given by its two-character tag name.
bam_aux2‹C› converts the pointer returned by the previous function to the actual type of the optional tag; here, ‹C› is a type identifier: i for an integer, Z for a zero-terminated character string, etc. Refer to the documentation and the htslib/sam.h file for a complete list.

  Does it keep the order of the optional fields, or is this information unavailable?

As long as you don’t modify optional fields, the order is kept, but unspecified. Do not rely on the order!

  When a user runs htslib or even samtools view, how does the BAM row "end"

BAM records store a header that contains the record size.


Convert or extract Seurat object as the input of FateID I am currently using Seurat to do the scRNA-seq clustering analysis. After this, I am planning to use FateID to do the downstream analysis, could someone teach me how to convert or extract a Seurat object as the input for FateID?

Thanks in advance.

Seurat objects have various slots where different types and forms of data are stored. For FateID you would need an expression matrix (with cells as columns and genes as rows) and this can be extracted via:

seurat_object@assays$RNA@data

Here RNA is the assay (you can have multiple assays per sample, RNA, protein, ...) and data is normalized data. You can extract counts or scaled data instead.

If you are using Seurat v2, drop the @assays$RNA part.


Add segregating sites to branches of a tree I'm trying to figure out how I would plot a tree with number of segregating sites on display on the branches
I'm using acctran from phangorn to plot a parsimony tree (using phydat object from fasta files)
I have been trying to use seg.sites from the ape package to somehow merge these objects together but there has to be some easy way to plot a tree with the information that it was made from

The segregating sites in parsimony are formally called steps. Any number of programs will do this but I don't really use R, particularly for phylogeny, You can of course simply look at the raw treefile, all the data is there (described in the last paragraph). If look at the options for genetic distance, in parsimony this will be an integer and again the technical name is steps and is the segreating sites you want.
Python's ETE3 package will work (but you'd need to pipe it through Biopython), as well as specific standalone packages Figtree, the steps or 'genetic distance' on the branches on the branches option is in the check box. The heavy weight program for tracking step associated phylogeny using parsimony is Mesquite and its famous grandparent is called MacClade.
I suspect its a nomenclature thing, the options probably talk about genetic distance because few people do parsimony phylogeny now. Thats the option you want  and you will know because it will be >>1, whilst distance and model-based methods will be between 0 to 1. From a tree viewing it doesn't matter if the genetic distance is 0 to 1 or 1 and beyond because its all proportional. Its just an option thing you are missing. Thus if there is an option, e.g. to include bootstrap support on the tree, there'll be an option for genetic distance.
Just make sure you don't confuse the number of steps per branch with a consensus statistic, because there're completely different but might give similar looking outputs. Again its easily checked against the raw treefile because the "steps" will be immediately adajacent to the taxa name for the "twigs" (terminal branches) in the tree.


How to find which scales are available in SeqUtil of Biopython? I am trying to analyse the hydrophobicity of a sequence using BioPython's "SeqUtil".

analysed_seq = ProteinAnalysis("AKVLGSLEALEPVHYEEKNWCE")
analysed_seq.protein_scale(Scale, WindowSize, Edge)

Known parameters

WindowSize refers to the total number of amino acids to be considered either side of the residue (an int variable).
Edge refers to the weighting of the scale further from the centre of the window (0-1). `

Unknown parameter

Scalerefers to which scale being used. I am interested in how to give theScale` variable a valid variable.

The Biopython docs simply point out that there are many scales to choose from and the source code doesn't seem to give any clues about scale options.

Where can I find a list of the scale are available and how to call them?

Use this line of code:

print(analysed_seq.protein_scale(analysed_seq.get_amino_acids_percent(), window=9, edge=0.4))


scanBam from Rsamtools is not importing one of my reads into R I have this read in my BAM file. It maps on chromosome 1.

I open this BAM file in IGV, and I can see the alignment on chromosome 1.

But when I open this file in R with Rsamtools:

bamContigsCel <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref     /bristolAssemblySorted.bam', param = Rsamtools::ScanBamParam(what = Rsamtools::scanBamWhat(), flag = Rsamtools::scanBamFlag(isMinusStrand = FALSE), tag = bamTags))[[1]]

I then check if the read maps to chromosome 1 in my R object but I cannot find it.

bamContigsCel$rname[bamContigsCel$qname == '000000F|arrow']
[1] II II IV
Levels: I II III IV MtDNA V X

But if I look at the BAM file, it's there. Why isn't Rsamtools importing my read into R?

$ samtools view bristolAssemblySorted.bam | grep -n '000000F' | head  -c80
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   0 0   GAAGCTGTCTAAACTTTGGC

Cross-posted on biostars

Note the flag; that read is mapped in a reverse-complemented manner, so isMinusStrand = FALSE is filtering it out. I tested this by making a BAM file with only that read:

@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

Then in R:

> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)


Error in FindIntegrationAnchors- Seurat package I am working on integrating a labelled single cell RNA seq cell atlas with an unlabelled one. I am wondering how do I determine the value of max.features when integrating two single cell RNA seq datasets?
Below is the code I am using:
epithelial = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/epithelial/filtered_feature_bc_matrix/",gene.column=1,cell.column=1)
habermann = Read10X(data.dir="/projects/b1025/sdi0596/Covid19_Single_Cell_RNASeq_project/filtered_matrices_count/Habermann_Single_Cell_Data/",gene.column=1,cell.column=1)

data.list=c(epithelial,habermann)

normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Perform integration
immune.achors= FindIntegrationAnchors(object.list = ifnb.list,max.features = 200,verbose = TRUE)

Error in FindIntegrationAnchors(object.list = ifnb.list, max.features = 200, :
Max dimension too large: objects 1, 2 contain fewer than 30 cells.
Please specify a maximum dimensions that is less than the number of cells in any object (4).

Regarding the dimension of these objects, the epithelial one is 33538 x 49109 and the habermann is 33694 x 220213.
Any insights would be highly helpful.
Thank you.

There's a few problems with your code, first, when you do Read10X() it returns you a sparse matrix, and you need to put this into a Seurat object with meta data, before doing the integration.
So for example, i create some example data that is similar to your output from Read10X() :
library(Matrix)
epithelial = Matrix(matrix(rnbinom(100000,mu=200,size=1),ncol=50),sparse=TRUE)
rownames(epithelial) = paste0("gene",1:2000)
colnames(epithelial) = paste0("e",1:50)

habermann = epithelial + rnbinom(100000,mu=20,size=1)
colnames(habermann) = paste0("h",1:50)

Then convert to a seurat object:
data.list=c(epithelial = CreateSeuratObject(epithelial),
            habermann = CreateSeuratObject(habermann))

Then find variable features:
data.list <- lapply(data.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

Should look like this:
$epithelial
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

$habermann
An object of class Seurat
2000 features across 50 samples within 1 assay
Active assay: RNA (2000 features, 2000 variable features)

Now perform integration, below I have to reduce k.filter because I have very little cells in this example. In your case, you can simply use the default settings.
features <- SelectIntegrationFeatures(object.list = data.list)
anchors= FindIntegrationAnchors( data.list,max.features = 200,
k.filter=50,k.anchor = 3,verbose = TRUE)


Mapping statistics from bam file using bbtools and sambamba Below are the statistics for RNA-seq mapped and unmapped paired-end reads to rice genome using reformat.sh from bbtools on bam files. It gives 77% mapped and 5% unmapped, what about the remaining 18% reads? How I can get information for all reads?
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.mapped.bam, mappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.mapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 52020538 reads (77.10%) 7208830835 bases (77.41%)

Time:                         171.212 seconds.
Reads Processed:      67471k 394.08k reads/sec
Bases Processed:       9312m 54.39m bases/sec

For unmapped read;
Command;
reformat.sh in=Leaf_T1_F_R10_S1_L001.bam out=Leaf_T1_F_R10_S1_L001.unmapped.bam, unmappedonly

Unspecified format for output Leaf_T1_F_R10_S1_L001.unmapped.bam,; defaulting to fastq.
Found sambamba.
Input is being processed as unpaired
Input:                   67471075 reads           9312399845 bases
Output:                 3435326 reads (5.09%) 465456267 bases (5.00%)

Time:                         142.978 seconds.
Reads Processed:      67471k 471.90k reads/sec
Bases Processed:       9312m 65.13m bases/sec

With sambamba; Can I say 94.91% reads are mapped and 5.09% are unmapped?
67471075 + 0 in total (QC-passed reads + QC-failed reads)
12015211 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
64035749 + 0 mapped (94.91%:N/A)
55455864 + 0 paired in sequencing
27727932 + 0 read1
27727932 + 0 read2
51006046 + 0 properly paired (91.98%:N/A)
51419222 + 0 with itself and mate mapped
601316 + 0 singletons (1.08%:N/A)
263816 + 0 with mate mapped to a different chr
237111 + 0 with mate mapped to a different chr (mapQ>=5)

Try samtools flagstat or sambamba flagstat to get similar type of information about your bam file. It may be easier to interpret. Both of these tools can be installed using conda if needed.
EDIT:
Yes, the results of sambamba flagstat you showed in your edited post indicate exactly what you said: 94.91% mapped and 100% - 94.91% = 5.09% unmapped:
64035749 + 0 mapped (94.91%:N/A)


Getting exact mismatch data from multiple sequence alignment TL;DR Can anyone point me in the direction of a tool that can be used to input different viral substrain reference sequences (fasta) as well as sequence data (fastq/bam/vcf) and find the most likely substrain of the sequenced data?
I am doing multiple sequence alignment of different viral subtype references and would like to get the exact positions and base each time there is a mismatch. I have tried using different tools. For example Clustal Omega can almost give me what I seek, but I have trouble getting the exact position and nucleotide for the mismatches (the full thing is ~8000bp long):
ACTACAATAATCCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAATTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG

*********** *********** ****** *****************************

E.g I would like something like:
Pos 12:
C
T
T
T
T
T
T
T
T
T
T

For each position with a mismatch.
A problem is that the sequences are slightly different lengths so something like a variant call format on the "main" strain might work, but I'm not sure how I can do this on multiple files of this length.
The end goal is actually to be able to call which substrain is present from NGS data. I imagine using variant calling on the "Main" strain, but getting from there to calling which substrain is in the data is a bit of a long way for me. If someone have any suggestions on tools for this, please enlighten me :)
Thank you very much!

I have partly answered this for know! For anyone interested i use kalign to align all 10 strains to a fasta format, then I use snp-sites to output a vcf file that can show each mismatch in accordance to the first strain in the fasta. New problem is that the alignment "stretches" the length of this first strain, because some of the other strains are a bit longer. This means that the position of variants are with reference to the longest strain of the 10 and not the strain I want. It can work for now, but any solution to this, would be helpful.


Is there a simple command for outputting a tab delimited columns? I am working on a fasta file and am writing my command in nano within command-line and executing using python, also within a command line.
My objective is to get my command to provide me with a tab delimited file with three columns: first column should contain my sequence name, second column should provide me with my sequence length, and the third column should show the sequence itself.
I have written the following command so far within nano:
For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq
A06842X 45 GORHQRIHQEERHEUWOHRPPTRWFAWWEAKJNFWEJQEFQEPRT

Jb4329PC 21 LDHKQAJTYOWEPLKFREEQW

I have no idea what do you want to do with fasta files reformatted like this, but I think something like this should do the job (untested):
from Bio import SeqIO
import sys

for olig_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  name = olig_fasta.name
  seq = str(olig_fasta.seq)
  seq_len = str(len(olig_fasta.seq))

  print("\t".join([name, seq_len, seq]))


Matlab scatter diagram I have two columns of data, x and y, how can I make a ‘thermodynamic diagram’ with matlab? The z axis is the density value of dots, like below figure.

The x-axis is angle, y-axis is distance.
I want to construct a heatmap like this figure:

A reverse image search on google leads to the exact code used to produce the second image:

[Xcoarse, Ycoarse] = meshgrid([0 1 2 3], [0 1 2 3]);
[Xfine, Yfine] = meshgrid(linspace(0,3,3000), linspace(0,3,3000));
DataCoarse = [ 1 2 4 1; ...
           6 3 5 2; ...
           4 2 1 5; ...
           5 4 2 3];
DataBicubicFine  = interp2(Xcoarse, Ycoarse, DataCoarse, Xfine, Yfine, 'bicubic');
figure
surf(Xfine, Yfine, DataBicubicFine); shading flat; colormap(jet); view(0, 90)
 hold on; plot3(Xcoarse, Ycoarse, 10*ones(size(Xcoarse)), 'k.', 'MarkerSize', 20)
colorbar
print -dpng -r300 BicubicInterpolationExample.png


Purpose of correlating docking score with catalytic efficiency(Kcat/km) What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?
I saw some papers are doing this, is it to validate the docking protocol?
Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

What is the purpose of correlating docking score with catalytic efficiency (kcat/km)?

This is specific to the enzymology of paper you are citing. If the solubility of your substrate is lower than your Michaelis constant (K_M) you cannot fit to a Michaelis-Menten curve and can only get the slope of the initial straight part of the curve, which is the catalytic efficiency (k_cat/K_M). Ideally K_M would be better. But in generally catalytic efficiency ought to correlate with well calculated ddGs of ligand binding.

  I saw some papers are doing this, is it to validate the docking protocol?

Totally. Docking finds were ligands bind. Although the energetics don't match too well the experimental kinetics as of 2019.
The catalytic efficiency (k_cat/K_M) is the ratio of the turnover (k_cat) over Michaelis constant (K_M). The Michealis constant is in turn the ratio of on-rate over off-rate plus turnover. The ΔΔGibbs of the ligand ought to correlate with the ratio of the on-rate/off-rate, while the turnover is unrelated.
The turnover should correlate with the difference in energy between the substrate bound form and the transition state. This is the realm of theozyme design and requires correctly protonated ligands and protein, but when it is you get a nice transition state hill (e.g. personal example).

There are a few "ought to" in the above. Docking and force-field energy scores are not precise. The paper you refer to is doing very advance docking (akin to dynamic undocking) and, not accidentally, the ligands are very flexible and rather soluble.

  Lets say the docking protocol is correct, the docking score will positively correlate with kcat/km?

A positive correlation coefficient just means that generally both the catalytic efficiency and docking scores increases or decrease together. If a metric "worked backwards", such as ΔΔGibbs (more negative is better), you would want a negative correlation. What you want is a strong correlation, with the sign a result of the nature of your metrics.


How can we distinguish between true zero and dropout-zero counts in single-cell RNA-seq? In single-cell RNA-seq data we have an inflated number of 0 (or near-zero) counts due to low mRNA capture rate and other inefficiencies.

How can we decide which genes are 0 due to gene dropout (lack of measurement sensitivity), and which are genuinely not expressed in the cell?

Deeper sequencing does not solve this problem as shown on the below saturation curve of 10x Chromium data:

Also see Hicks et al. (2017) for a discussion of the problem:

  Zero can arise in two ways:

  the gene was not expressing any RNA (referred to as structural zeros) or
  the RNA in the cell was not detected due to limitations of current experimental protocols (referred to as dropouts)

Actually this is one of the main problems you have when analyzing scRNA-seq data, and there is no established method for dealing with this. Different (dedicated) algorithms deal with it in different ways, but mostly you rely on how good the error modelling of your software is (a great read is the review by Wagner, Regev & Yosef, esp. the section on "False negatives and overamplification"). There are a couple of options:

You can impute values, i.e. fill in the gaps on technical zeros. CIDR and scImpute do it directly. MAGIC and ZIFA project cells into a lower-dimensional space and use their similarity there to decide how to fill in the blanks.
Some people straight up exclude genes that are expressed in very low numbers. I can't give you citations off the top of my head, but many trajectory inference algorithms like monocle2 and SLICER have heuristics to choose informative genes for their analysis.
If the method you use for analysis doesn't model gene expression explicitly but uses some other distance method to quantify similarity between cells (like cosine distance, euclidean distance, correlation), then the noise introduced by dropout can be covered by the signal of genes that are highly expressed. Note that this is dangerous, as genes that are highly expressed are not necessarily informative.
ERCC spike ins can help you reduce technical noise, but I am not familiar with the Chromium protocol so maybe it doesn't apply there (?)

since we are speaking about noise, you might consider using a protocol with unique molecular identifiers. They remove the amplification errors almost completely, at least for the transcripts that you capture...

EDIT: Also, I would highly recommend using something more advanced than PCA to do the analysis. Software like the above-mentioned Monocle or destiny is easy to operate and increases the power of your analysis considerably.


Identifying proteins in a concatenated protein sequences I have a 120-concatenated protein sequence in a fasta file. I would like to split by proteins.
How can I identify each protein?

The only way I can think is to break the protein apart on methionine residues, blast each amino acid locus and concatenate any locus that does not score a query protein position of '1'. The theory is simple every protein starts with a methionine residue, however there are plenty of M residues within the protein that don't represent the beginning of the gene.

Breaking the protein apart on M would be easy, but then automating and parsing the blast would be fiddly. Obvioulsy you will need to code it and if its a one off would it be worth it?


Primer-Blast Custom Database [error] I'm using Primer-Blast to compare a PCR template against a custom database of reference sequences, using Primer-Blast's Custom database functionality.

However whenever I try to upload a FASTA file as the custom database and then run Primer-Blast, I get the following error:

Exception error: Fasta Reader: sequence id ends with 80 valid nucleotide characters. Was the sequence accidentally placed in the definition line?

No matter how I try to format my FASTA file, I get this error.

Does anyone know how to properly format a FASTA file to have it accepted as a custom database of reference sequences to compare against the PCR template? Or is this a bug in the software?

Here's an example FASTA file, to upload both for the template and the custom database.

>test1
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTGC

I just copied your FASTA sequence into a file (on Mac OS) and uploaded it to Primer BLAST, the search did not give an error.

If the error is caused by the line endings, maybe you can fix your downloaded files with a tool like dos2unix.

Update:

On a whim, I tried adding a description to the FASTA header:

>test1 description
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTG

This time, there was no error.

Does that work for you as well?


Alternative to enrichR for enrichment analysis? I have been using EdgeR from the Bioconductor package to perform differential expression analysis on bulk RNASeq data sets.
This analysis gives me an output dataframe that has a list of all the genes, ranked by the their adjP values. This dataframe also contains the t-statistic and lfc values for each gene as well.
I have recently been using enrichR to perform gene set enrichment analysis on all the differentially expressed genes to determine what biological, molecular, transcriptional, etc pathways these genes are involved in. I really like the enrichR package for a few key reasons:
1.) It has access to many different databases that I can call on with a simple vector of their names. Most of my data is murine RNASEq, so for those analyses, I can simply select the murine databases:
enrichR_dbs <- c("WikiPathways_2019_Mouse", "KEGG_2019_Mouse", "GO_Biological_Process_2018", "GO_Cellular_Component_2018", "GO_Molecular_Function_2018", "GWAS_Catalog_2019", "TRRUST_Transcription_Factors_2019")

2.) It uses a single function, enrichR, to perform the enrichment. The arguments for that function are the list of gene symbols and the list of databases.
enrichr(as.character(Gene_Dataframe$SYMBOL[1:500]), enrichR_dbs)

Anyway, I am looking for another enrichment analysis because I want to be able to see if I get similar pathways for my genes when using another method/algorithm.
Unfortunately, I can't seem to find another enrichment analysis( competitive or self-contained) that is both as user-friendly and that gives me easy access to murine databases (without requiring to convert my genes through changes in capitalization, like you need to with GSEA) as enrichR.
Does anyone have recommendations about enrichment analysis R packages that don't require too many more inputs than the dataframe I generated through my EdgeR analysis?
P.S. I'm ok with downloading databases if I need to. Though inclusion of a PDF with a tutorial would be fantastic, if possible.

A couple of years ago I wrote a blog post about GSEA methods in Bioconductor, in addition to these I would recommend clusterProfiler or ReactomePA.  In most of them you need to provide the gene set/pathway database.
If you want to do gene ontology enrichment then I recomment the topGO package.
In most, if not all, of these methods you just need an ordered list and a gene set source.


gene-gene correlation from two different Tissues I am stuck on how to do correlation for two independent data sets with common row and column names.
A and B are datasets that contain as many rows as genes and as many columns as samples.
The rows in A and B represent a common set of genes but measured in two different tissues.
The columns represent measurements in the same 5 samples in both A and B.
I want to do a correlation between the set of genes in A and B. This is to see if the same genes in both tissues are correlated or not.
Since the matrix would be big in my actual data, I only want to retain a correlation coefficient higher than 0.5.
Here I simulate the data set.
set.seed(1)
A <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(A) <- paste0("G_", 1:100)
colnames(A) <- paste0("M_", 1:5)

set.seed(42)
B <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(B) <- paste0("G_", 1:100)
colnames(B) <- paste0("I_", 1:5)

Thank you!


run a tool with python>=3.6 on galaxy-project instance that runs on python2.7 I am developing a wrapper for a bioinformatics tool writen in python3.

The native Galaxy-project runs on python2.7.

How can I the tool that requires python3 in a galaxy-project instance?

Each Galaxy tool uses its own (typically conda-based) environment, which will normally not use the same python as Galaxy itself (note that Galaxy can run under python 3). An example of this is deepTools, which is python 3 only but will happily run under a python 2 Galaxy since conda is being used to create a python 3 environment for it. It's best to distribute your tool via bioconda and then Galaxy will be easily able to install it.


Count files using htseq-count? I am doing RNA-Seq data analysis using data available on SRA (SRR6047326). Samples are from the pig and they are paired. So, I did alignment using HISAT2.I used pig reference genome available on ENSEMBL. Then, I sorted .sam output files obtained from HISAT2. Now I am trying to get the count file using sorted sam file and it is not working. I am using annotated GTF file for pig which I downloaded from Ensembl.

I get count files using following commands, but all the counts are 0.

htseq-counts -i no myfile.sam genes.gtf>output.count

What is wrong here ? Are there any alternatives to htseq-count?

I m not sure about htseq-counts but I would suggest go for featureCounts, I will update my answer and share the script/command. Give it a try

the command i use

featureCounts -T 50 -C -p -t exon -g gene_id -a /run/media/data2/gencode.v21.annotation.gtf  -o PRIMARY_CELL.txt $(ls *.bam | sort )


Performing PCA for the samples and for the genes I have 10 samples from a RNAseq experiment (5 control, 5 disease), I have performed a cluster analysis for the samples and for the genes (4000 genes aprox) to see how they cluster (to see which samples are similar and which genes have similar expressions).
I was wondering if in terms of statistics it makes any sense to perform a PCA of the samples instead of the genes?
thank you in advance!

First of all, I am not a statistician but have been performing RNA-seq analysis for a while, here is my take:

I think performing PCA on the samples makes sense in terms of mathematics. What I mean is you can compute eigenvectors and subsequently principal components either on your genes or on your samples.

I don't think performing PCA on the samples makes sense in terms of statistics in RNA-seq context (and probably in any other context): Your genes (features) describe your samples and not the other way around. So the computed PCs should describe/explain your samples and not the other way around. You would not perform PCA on your samples (individual cells in scRNA-seq) even in the case of scRNA-seq, where number of samples sometimes exceeds the number of genes, there PCA is performed on genes and the resulting few PCs are used for downstream analysis, usually clustering and visualization, both used to infer about samples (cells).


Database for proteome-wide predictions of protein structures Accuracies of protein structure predictors have improved quite a lot in recent years. Algorithms such as Rosetta have gotten robust enough to predict structures of large number of proteins. However, I could't find any initiative to make a database for proteome-wide predictions of protein structures or may be integrate the predicted structures in existing databases (PDB or Uniprot).

If any such exists, please let me know.

There are three initiatives I know of to have a go at this:

PConsFam, which collects data from this paper.
The Baker group's metagenomic study which you mention in the question.
The recent DMPfold work from our lab. This compares its results to the above two studies and discusses the effect on various model organisms.

The second and third of these don't provide a database per se, but the models are available for download and are linked to the Pfam database.


hgsql not found when running a script despite being accessible on command-line I am working on macOS. I am following the steps described here to lift over an annotation from the reference genome of C. elegans to another assembly of C. elegans.

I have run this code:

time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

However, the do.log file contains:

zsh:1: command not found: hgsql
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Sun Nov 11 15:38:41 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/job.csh
# chmod a+x /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
# ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh
cd /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat
rm -rf tParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 58018742 0 /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.chrom.sizes 2000 -lstDir=tParts
rm -rf qParts
/Users/username/Documents/projects/main/output/software/liftover/data/scripts/partitionSequence.pl 10000000 0 ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes 1000 -lstDir=qParts
Can't open seq size file ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.chrom.sizes: No such file or directory
Command failed:
ssh -x -o 'StrictHostKeyChecking = no' -o 'BatchMode = yes' localhost nice /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/run.blat/doAlign.csh

I think the biggest problem is that the hgsql command is not found.

If I run hgsql on the command line, I find it even though I get the following error:

dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
  Referenced from: /usr/local/bin/hgsql
  Reason: image not found
[1]    34059 abort      hgsql

Location of the hgsql installation:

which hgsql
/Users/username/Documents/projects/main/output/software/liftover/data/bin/hgsql

Can anyone help me with using hgsql, please?

It is possible that the shell being used in the process is different from your standard command line shell. I am not familiar with the tool, and it's hard to say more without knowing your OSX version, but it looks like it is using zsh. I think that zsh only fairly recently became the standard OSX shell. I don't know the tool or the shell very well, but it looks like you might be able to explicitly set up hgsql in your PATH in your .zshrc.

It's also possible that the liftover script was not built to run via zsh. I notice that it is trying to run csh scripts as part of the workflow. Then again the localhost may insulate it.

So I think that my first step would be to sort out which shell I am using, whether it works with the liftover pipeline, and how it's running on the localhost. I notice that the ssh call running the csh script on the localhost is the step that actually fails. What shell is executing the csh script? Does that shell know about hgsql? It could be that the way that the error bubbles up it is hiding which shell it is actually using.

I'd recommend trying that ssh call in isolation from the rest of the pipeline, and diagnosing the error that way.

Hope that helps.


Error Viewing file on HiGlass I am implementing hi-glass python on Jupyter, and upon running the test data code

from higlass.client import View, Track
from higlass.tilesets import cooler

ts1 = cooler('../data/Dixon2012-J1-NcoI-R1-filtered.100kb.multires.cool')
tr1 = Track('heatmap', tileset=ts1)
view1 = View([tr1])
display, server, viewconf = higlass.display([view1])

display

I get the error message stating:

Error Retrieving tileset info: http://localhost:60501/spi/v1

Please advise on the correcting. The code above has been copied from the HiGlass webpage.

I'm one of the authors of HiGlass. Sorry you're getting that error. I just tried it out in a fresh environment and it worked for me. Just one thing is missing in the docs which is that you need to do import higlass.

Here's what I did to set up my environment:

conda create --name test-bi python=3.6
source activate test-bi

The rest is from the docs:

pip install jupyter higlass-python
jupyter nbextension install --py --sys-prefix --symlink higlass
jupyter nbextension enable --py --sys-prefix higlass
jupyter notebook

It looks like either you or somebody else added a bug report here: https://github.com/higlass/higlass-python/issues/29

I'll respond and continue the conversation there.

Also, there's a repository full of examples from the HiGlass SciPy presentation here: https://github.com/higlass/scipy19


sort a fasta file containing the Oxford Nanopore Technologies (ONT) header by sequencing start_time ascending I have basecalled ONT reads and converted them to multifasta. The multifasta contains the original ONT headers in this format:

3ebb56cd-3671-4abd-b1ac-0c759bt068d0 runid=eb6214851489c8e00eb0dcbd00d737f7dddxxxx read=1078 ch=422 start_time=2017-11-21T20:35:34Z

By checking the sequencing output as fasta, the sequences are not ordered by start_time.

>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC

I learned how to extract header IDs and time_start from ONT fastq files by using this answer.

I'd like to have a similar approach but with fasta.

Yet, I would like to create a new fasta file where the sequences are written in an ascending order as they were sequenced (start_time) as follows:

>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG

Edit: After doing some reading, if it can be assumed that all the timestamps are ISO8601 and of the same representation of ISO861 (there are a few), then they can be sorted lexicographically, directly, without converting to epoch-time integers. So I have edited this answer to add an answer without conversion steps:

Here's one approach that uses the timestamps directly:

(1) Extract the start times from the headers and sort them to another file, in ascending order:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa | sort > starts_iso8601.sorted.txt

(2) Use that file of ordered ISO8601-formatted starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Quite a bit cleaner and faster.

You might separate the timestamp extraction and sort steps into two commands and inspect the intermediate result, just to ensure that you have valid strings. Or use the approach below to convert the timestamps, which will catch bad input.

Here's a second way that may add unnecessary work. As an advantage, however, it adds some error checking by converting the ISO8601-formatted timestamp.

If the string you pull out of the header is not a valid timestamp, the first conversion step (ISO8601-to-epoch) will fail and you will then be able to catch invalid input pretty quickly, which can add confidence that your output is correct.

(1) Extract the start times from the headers:

$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa > starts_iso8601.txt

(2) Convert the ISO8601-formatted times to time-since-epoch times:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"$line" +%s; done < starts_iso8601.txt > starts_epoch.txt

(3) Sort the epoch times — which are just integers — to another file, in ascending order:

$ sort -n starts_epoch.txt > starts_epoch.sorted.txt

(4) Convert sorted epoch starts back to ISO8601 format:

$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"@$line" -u +"%Y-%m-%dT%H:%M:%SZ"; done < starts_epoch.sorted.txt > starts_iso8601.sorted.txt

(5) Use that file of ordered starts to write out FASTA records in your sort order:

$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

Note: Use gawk, gdate, gsort in place of awk, date and sort if you're using OS X, Homebrew, and GNU core utilities.

I'm showing the steps separately to explain the general procedure.

Many or all of these steps can be glued together in a Unix pipeline to save the time, memory and disk space costs of writing intermediate files.

Also, I'm not doing much error checking here. You might look at cases where headers do not contain start_time, or times are not all in ISO8601, or other cases not anticipated here.


Is it possible to merge scRNAseq data from experiments with different design? I have 4 different single-cell RNAseq experiments, each one representing a different sample of cell type population. I'd like to merge them to a single dataset. However, different cell types are enriched in each sample, and different samples can represent also different conditions (control and disease).

If I run a PCA on the merged dataset (without any processing of data) I find that the different sample's cells cluster separately. This could be due to biological and technical (i.e. batch effects) variation.

Is it possible to merge these different samples while retaining a good amount of information or should they be analyzed separately?

It’s right there on the cellranger manual:

#aggregate results of counts for separate samples
cellranger aggr
#analyse the combined results
cellranger reanalyze

Note: I’m not sure this applies to the above case where batch effects could be an issue. However, this is how the analysis is performed for 10X data from multiple runs in principle.

cellranger aggr

cellranger reanalyze


WGCNA error codes in Network Construction & moduel detection I am using WGCNA and code has worked with no problems with some of my datasets however with most recent attempt when completing gene network & identification of modules. I get following error code...

net__30= blockwiseModules(t5, power=6,
                              TOMType = "signed", minModuleSize=30,
                              reassignThreshold =0, mergeCutHeight = 0.25,
                              numericLabels = TRUE, pamRespectsDendro = FALSE,
                              saveTOMs = TRUE,
                              saveTOMFileBase = "PBCreponseTOMB",
                              verbose= 3)

I come up with error code saying:

            Number of items to replace is not a multiple of replacement length.

Does anyone please know what this indicts to? I have looked on all the WGCNA tutorial pages but cannot find this error message anywhere.

Thank you kindly

Looks like a bug in the code which I will try to hunt down and fix. In the meantime, I would suggest playing with maxBlockSize argument to blockwiseModules. Try increasing it as much as your available RAM allows (see the paragraph "A second word of caution concerning block size" on page 6 of WGCNA tutorial I, section 2c, for some guidance on how to set maxBlockSize). Hope this helps.


Obtaining Whole Genetic Sequence As an end-user of my own data, specifically raw DNA sequence (WGS, whole genome sequencing). How or where do one obtain such DNA data so that a biology hobbyist can perform bioinformatics on this?
What keywords is best describe to the search engines when finding these online WGS websites).
It would be nice to have a USB stick mailed to me after such DNA HomeKit.  What file format is most commonly offered? (VCF?)
Am I looking too far ahead into the future?

"DNA data" comes in several forms dependent on the technology used to produce them.
Companies like 23andMe are using SNP chips and those are available for an only rather a limited number of species. To be completely honest, I don't even know what would be the raw data they would provide, but I would suspect it could be a vcf file.
If the chosen technology would be DNA sequencing. DNA sequencing is provided by a few dozen of sequencing centres around the world, some of them associated with universities. DNA sequencing generates files that contain individual sequencing reads with assigned quality scores ("how much the sequencer trusts each base"), the typical formats are fastq (Illumina) or .bam (long read technologies). Usually, you would also need to find a suitable reference in one of the public databases, for instance, the NCBI portal for human reference. Then map your reads and do variant calling. Of course, then it depends on what would you be interested in.
If you don't care for your own DNA and you are more interested in life around you, you can check ncbi for all sorts of genomic data. If you would be interested to understand how the raw sequencing data work, ebi is a good place to look. Also, in every single paper generating/using sequencing data you there should be mention where in public databases you can find them.
Finally, let me warn you. Some of the genomic analyses are quite computationally heavy (not all though). If you don't have a specific question in mind, it would be wise to start with bacterial sequencing.
Hope you will have fun with DNA data!


Get Dockerfile of a bioconda package Bioconda packages have associated containers images. E.g. bioconductor-csaw has the image quay.io/biocontainers/bioconductor-csaw.
I would like to customize the image, in particular, I would like to add rstudio.
It is possible to get the Dockerfile associated to the image of a bioconda package?
How do you suggest to proceed to customize the image?

The Dockerfile is generated by galaxy-utils, so I'm not entirely sure what it looks like. However you should be able to generate such a "mulled" container as documented here. In short, something of the form:
mulled-build build 'bioconductor-csaw=someversion--somebuildnumber' 'someotherpackage=someversion--somebuildnumber' ...

That probably requires setting the CONDA_IMAGE environment variable. Bioconda currently sets that to quay.io/dpryan79/mulled_container:latest.


Comparing factor levels in deseq2 I have sample information (sampleinfo) in this format:

SampleID    condition   type
MP_1    A_1 Paired-End
MP_2    A_1 Paired-End
MP_3    A_1 Paired-End
MN_1    A_2 Paired-End
MN_2    A_2 Paired-End
MN_3    A_2 Paired-End
WP_1    B_1 Paired-End
WP_2    B_1 Paired-End
WP_3    B_1 Paired-End
WN_1    B_2 Paired-End
WN_2    B_2 Paired-End
WN_3    B_2 Paired-End

I want to compare A_1 with A_2, B_1 with B_2 ans also want to compare (A_1+A_2) with (B_1+B_2).

This is how I used DESeqDataSetFromMatrix

ddsMat <- DESeqDataSetFromMatrix(countData = seqdata[,rownames(sampleinfo)], colData =  sampleinfo, design = ~ condition)
ddsMat <- DESeq(ddsMat)

res.05_FC0<-results(ddsMat,contrast = c("condition","A_1","A_2"),alpha=0.05,lfcThreshold = 0)
res.05_FC0<-results(ddsMat,contrast = c("condition","B_1","B_2"),alpha=0.05,lfcThreshold = 0)

How can I compare (A_1+A_2) with (B_1+B_2) using the contrast argument in results function?

Of course I can make a new sampleinfo table, but I want to use the same sampleinfo as above. I tried to read about numeric contrast, but I am unable to use it correctly.

You can use contrast = c(0.5, 0.5, -0.5, -0.5), which means "the average of the first two coefficients minus the average of the second two.


What unit (TPM, FPKM/RPKM, or other) to use when working across samples I have raw gene read counts and would like to perform an analysis across multiple samples. I've found conflicting info online on how this should be done. One commonality however is that FPKM/RPKM should not be used. On the other hand, there seems (to me at least) to be conflict about whether TPM is suitable. For example, this StatQuest video seems to suggest that TPM is suitable for analysis across samples. One person on this Biostar post (among other resources I've looked at) suggested that none of the units (FPKM, RPKM, TPM) are suitable for cross sample comparisons.
I am quite confused. Can someone clear up the confusion?

Neither of RPKM, FPKM or TPM is a good choice. All these methods are similar, they only perform the operations in a slightly different order. Eventually they correct for gene length and sequencing depth. Correction for sequencing depth is necessary but not sufficient in many cases. Correction for length is typically not done for a differential analysis as it reduces the counts of many genes and therefore throws away statistical power. I will not describe this in detail since the StatQuest video series you link includes videos about DESeq2 and edgeR normalization procedures which extensively cover the normalization procedures and why these approaches are superior to RPKM/FPKM/TPM. In short: You have to correct for differences in library composition to compare between samples. Please watch these videos and read the manuals of the established RNA-seq analysis packages such as but not exclusively limma, edgeR and DESeq2.


Protein model quality assessment How do I perform protein model quality assessment?

I will obtain predicted structures from DMPfold, I-Tasser, and Rosetta, what should be the next approach to select a structure from among this set of models?

Note, I have been looking at some, appropriate methods such as MetaMQAP, QA-RecobineIt, but their servers are not working.

Have a read of the Estimation of model accuracy in CASP13 paper to get the latest rankings.

Generally you would run the model quality software on the models you have generated and choose the best. Also have a look at the structures in PyMol, you learn a lot that way. See, for instance, how much variety there is, and where the models vary. Is there a conserved core in many of the models? Or is there just a blob of "spaghetti" that doesn't look like an ordered protein?

If you have many models (e.g. from a Rosetta run) you might want to look into clustering the models first to reduce the number you have to work with.


factoextra: Error in if (xlab == "x") xlab <- "x value" : argument is of length zero I just followed the tutorial trying to do clustering on my own data, but I can't visualize it because fviz_cluster function is giving me an error:

  Error in if (xlab == "x") xlab <- "x value" : argument is of length zero

Here is the code I am running:

library("fpc")
db <- fpc::dbscan(rtsne$Y, eps = 0.15, MinPts = 5)
# Plot DBSCAN results
library("factoextra")
fviz_cluster(db, data = rtsne$Y, stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
         geom = "point",palette = "jco", ggtheme = theme_classic())

How can I make fviz_cluster function working?

  Update

Filed a bug: https://github.com/kassambara/factoextra/issues/62

I actually found an answer just accidentally. It is very unfortunate that factoextra documentation does not explicitly say that data parameter should be of data.frame type. In my case rtsne$Y was a matrix. The error that factoextra is giving is very uninformative.

Here is the code that ultimately worked:

fviz_cluster(db,
----->  data = as.data.frame(rtsne$Y),
        stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
        geom = "point",palette = "jco", ggtheme = theme_classic())


Huge variability in blastn speed between same short sequence job I have been running a couple of blast searches and I am very suprised and somewhat frustrated by the huge variability in job completion time. The specifics of my blast search are as follows:

Running blast using webbrowser (chrome) via ncbi
19 nucleotide query
database: WGS of one particular sequencing project/bacterial strain (so database is VERY small, in total approx 4.5 MB divided over 60 contigs)
parameters are mostly untouched and are automatically optimized for short queries by blast

In one instance, I started the job more than an hour ago and it is still not done. Knowing from experience that this previously hasn't taken that much time, I started a new job in a different browser tab with exactly the same sequence, database and parameters, and now it finished within 20 seconds. What on earth is happening here?

First of all, if you're using an online service like NCBI's BLAST, then the time to completion will depend on the load the service is experiencing. Your job is ut in a queue and will be started when it's its turn. If the service is under heavy load, then its turn might come later than you expect. Then, if the service is under heavy load, even after launching your job, it might take longer to finish because it won't have enough resources to run quickly.

Since you say that the same job finished quickly on a new tab, it sounds like either the first job was last in the queue or there was just a bug that caused it not be launched or something like that. We can't really help with things that resolve themselves like that.


RefSeq transcript naming - transcript vs gene I have a dataset composed on numerous RefSeq transcripts with names of the form NM_###### - essentially standard annotated mRNA accession numbers.  I notice that numerous transcripts which have different sequences (different start and stop sites in a gene) all have the same NM_#### accession number.  The interpretation is that they are all from the same gene, but this makes no sense to me.  Why should different transcripts with different sequences have the same name?  I understand this identifies their gene of origin, but they are still different transcripts.

Here is an example from the .bed file I'm using.  Note they have different start and stop locations and even have different strand notations, but all have the same name.

chr     start        stop       name            score   strand
chrX    3441730     3443690     NM_001085545    0       -
chrX    3700232     3702192     NM_001085545    0       +
chrX    3750925     3752885     NM_001085545    0       -
chrX    31336368    31338328    NM_001085545    0       +

No, this is the same sequence. A single RefSeq accession will correspond to one sequence. However, nobody said that sequence needs to be unique in the genome. The specific sequence you are talking about, NM_001085545.1, is a mouse transcript which can be aligned to many different regions of the mouse's chromosome X.

This is your sequence:

>NM_001085545.1
agagaagtctaggcaggactggcagactacagcaatgcctcaaccaagaagcctctgcga
gagcttgttctaggaattttcatcaatggggcttttagtcagcagggtcttgagatgcag
ggattccagtctgttagagccacagccagaagccatagccggagccagctacattcctgg
cagtcgcaagcgaaaaagaaacagtttggaggagttggcaacaagttctaatgttcatgg
acctcaaaaccagggaatgtatccacatcaagttctcaactacatctactggaaaagggt
taagatctcatctaatgatgcttatcaaaacttatttttggacgggcatgatagcgacat
taaaatccgtgctctgggaaaaacatggtgtttacacaaagtatttttatgtcagtcagg
ctactttgctaacattctcaaaggtacttggagagaatcacaccatggtgttataaatct
gatcattaagaatgaggatattgatacccgatctctgcattttgtgtttggtgctttgta
cacggatgcggatttgtcaataacacctctggaagttcctcaagttttggcagcagcatg
cctgcttcgggtggatcaagtaattcagcagtgtgaaggaatcatgaaagaaactatcaa
caggaacactgtgtgctcctattatttggcagcagaaacctatagattaaaagctgtaaa
gacgagatgctttgaatggcttctttgcaatttgatggtacatccaagtgtggcacttta
caaggaagtagacttgaagttgatgtatcttctagcactgtcttctgacttactagtcat
gcaaaaggagattgatgtatataccacactaaaaatatggatgttcctttatcttaatcc
atgctggaacggaaccatgaaacagcttttacaacacgcaaacaactggctttccaccca
catggcatatgttgataacatcagttttcttgaaagtgaagaaggactaatatttcaacc
agtgtttaaaaagctgagatttcagcacattatctgtgacttgacttccacaactattct
tgaacaagatcgactaatacctatggcatggttgtcacccatttacaaacaacagtggtt
gactttgctgcgaacacaagaatatggggtaattggaccacaagttatcaatgaacaaga
acttgaagaatgcaccatgaggtgtggtacaatgatccccaaggatggaagatatacttg
gaagtggtcagttggacgacttggctttcctttacgtgtgacctttaccaggcagtgtgt
aattttaaggcaacggtgtcagaggtgtgatggttctgcttgccacaaccatatccgaaa
tgtcatattcagaataactttggtgtgttttgattccaacaaaagagtaactttcagaaa
gacaacaggttataaaatcctcacctttgaatataaggaggagcaaattgtaatgaaatt
ggatagtgatgttctaaccttccctatgtgtatattctgcaatttcctttttgtaaacct
aggaaatgcagaaaacaagtaatcttatcaattattagactgtaagtattttacaagtat
tgaatggttcatttaatatgaaggctgtgtaatgactacagtgaagaagaccaaaaacaa
attcaatattctttggagtctagtaccaccagaagtgactgcttccctccatcaccttat
ttccatgtatttctgaagaaaacctttacttccaacaaagtatttgtcaaatcaaataaa
cctaatcttctgatttaaatttaccatcgagatccaagcttggaatggaaataagagact
atatttgcaaggtaaaacatctcaatcaatttgttaagtagatttggccctagctaaaat
tgagacaaaataaaatcatgccattgttgaaactgttgtg

If you BLAT it against the mouse genome using the UCSC BLAT tool, you will see many hits, most of them almost perfect matches and some of them overlapping the regions given in your bed file:

Your bed file matches the following hits in the image above:

NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3441731   3443690   1960
NM_001085545.1  1960     1  1960  1960   100.0%  chrX   +     3700233   3702192   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3750926   3752885   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   +    31336369  31338328   1960

I don't know why your bed file has all of these various matches for NM_001085545.1, since you haven't told us what the bed file is supposed to be doing. A bed file is just a collection of regions. For all I know, that bed file is just "a few hits for sequence NM_001085545.1 in the genome". There is no reason to assume that each region in a bed file corresponds to a gene.

So there's nothing wrong with your file, you're just dealing with a non-unique sequence from the mouse genome. And you don't have multiple sequences with the same accession, you have one sequence, one accession, but multiple hits for that sequence in the genome.


How are positions visually indicated in the UCSC Genome Browser? I'm a little confused about how the position indices work in the browser. Based on the picture below, is base number 755440 an A or a G?

The general question is how do the indexing guides work? Is the base right under the | the one corresponding to the index or is it the one just before the |?

It's an A. The number refers to the position just before the |. This is clearer on different zoom levels, but I double checked by going to the very beginning of a chromosome:

As you can see in the image above, the | is always after the base it refers to. I also confirmed by checking a few random positions using another genome browser which places the positions unambiguously above each base.

For more details, see this UCSC blog post.


After KEGG and GO analysis, how to make tables+phylogenetic trees hope everyone is ok. I used Trinity to do a de novo transcriptome assembly, then blastp/blastx and then used Blast2GO software to do KEGG and GO analysis. So i got some txt files with header :

for GO

Sequence name Sequence desc. Sequence length Hit desc. Hit ACC E-Value Similarity Bit-Score Alignment length Positives

for KEGG

Seq Pathway Pathway ID #Enzs in Pathway Enzyme

All well so far. Then i was asked to do 3 tasks.

Tables with terpene synthases, (monoterpenes, sesquiterpenes, diterpenes and triterpenes) and corresponding phylogenetic trees including other known terpene synthases from all categories from other plant species.
Tables with Cytochrome P450s enzymes implicated in terpene biosynthesis (categorization per clades and clans) and corresponding phylogenetic trees including other known terpene CYPs from all clades and clans from other plant species.
From the graphs corresponding to the different KEGG and GO ontology analyses, export the data for the GO analysis for all three parameters (categories, names, blasts, total number of contigs and percentage) to an excel file to be able to transform them in other formats for future publication.

I have no idea how to proceed. Any guidance? Thanks

I'm going to answer 1+2 because 3 is hard to understand and maybe trivial?

Each of 1+2 seems to be have these sub-tasks:

identify the contigs in your assembly that have the genes of interest.
pull out the gene sequences from the contigs (you should have intermediate tblastx output somewhere with coordinates, if you don't you can just re-blast those contigs).
possibly translate those in silico to protein sequences.
tblastx (gene) or blastp (protein) those gene/protein sequences against whatever db is relevant.
pull sequences for the hits above some threshold for each gene from the db.
perform a multiple sequence alignment for each gene against the hits.
plug the MSA into whatever phylogenetic tool you like. I used phyml a lot back in the day, don't know what the kids like now.
make a visualization of the tree, probably?


Using shells other than bash As someone who's beginning to delve into bioinformatics, I'm noticing that like biology there are industry standards here, similar to Illumina in genomics and bowtie for alignment, many people use bash as shell.

Is using a shell besides bash going to cause issues for me?

Bioinformatics tools written in shell and other shell scripts generally specify the shell they want to use (via #!/bin/sh or e.g. #!/bin/bash if it matters), so won't be affected by your choice of user shell.

If you are writing significant shell scripts yourself, there are reasons to do it in a Bourne-style shell.  See Csh Programming Considered Harmful and other essays/polemics.

A Bourne-style shell is pretty much the industry standard, and if you choose a substantially different shell you'll have to translate some of the documentation of your bioinformatics tools.  It's not uncommon to have things like

  Set some variables pointing at reference data and add the script to your PATH to run it:

export FOO_REF=/path/to/stuff
export PATH=/path/to/foo-x.y:$PATH
foo blah blah

These will typically be shown in Bourne-shell syntax.  By using a different shell you have to translate the export commands to your local syntax, and especially PATH munging is somewhat shell-dependent.

If you're experienced in Unix, this will be only a minor niggle.  If you're a beginner, IMHO this will add a non-negligible amount of friction on top of all the other things you're learning.