Drawbacks of upper quartile normalization for scRNA-seq data I would like to use Upper Quartile normalization for scRNA-seq data defined as:

  The upperquartile (UQ) was proposed by (Bullard et al. 2010). Here each column is divided by the 75% quantile of the counts for each library. Often the calculated quantile is scaled by the median across cells to keep the absolute level of expression relatively consistent. A drawback to this method is that for low-depth scRNASeq experiments the large number of undetected genes may result in the 75% quantile being zero (or close to it). This limitation can be overcome by generalizing the idea and using a higher quantile (eg. the 99% quantile is the default in scater) or by excluding zeros prior to calculating the 75% quantile.

In particular, I'd use the 99% quantile to avoid normalizing by small number or zero. Is there any other drawback for using it on scRNA-seq data?

Bullard, James H, Elizabeth Purdom, Kasper D Hansen, and Sandrine Dudoit. 2010. “Evaluation of Statistical Methods for Normalization and Differential Expression in mRNA-Seq Experiments.” BMC Bioinformatics 11 (1). Springer Nature: 94. doi:10.1186/1471-2105-11-94.

My primary concern with using the top ~1% or so in upper quantile normalization is that it's going to be prone to the same robustness issues that RPKM/FPKMs run in to. That is, if for whatever technical reason you have to have a fair bit of variability in a couple very highly expressed genes (typically rRNAs, but one can imagine other genes) and the set of genes going into the UQ normalization is composed primarily of those then the normalization results will simply reflect irrelevant technical variability.

Of course, this all assumes that it's even appropriate to normalize scRNAseq data with a standard method like this. In some cases it is (e.g., when comparing control vs. treatment with a largely homogeneous cell population), but if you're dealing with very different cell types then I worry that the standard method is going to wash away things like differences in bulk RNA content between cell types. Presuming there are spike-ins, they should have the same normalization applied and then checked to ensure there's nothing greatly amiss. Having said that, if you're going to be doing tSNE and want to look for cell types, then odds are good you'll still be able to see most of them (after all, it's unlikely that RNA quantity is the only difference between them, though not removing that difference would make things easier).

The problem of global shifts in expression

The concept of global shifts in expression is not unique to scRNAseq, in fact the most common example of this is Myc up/down regulation in relation to cancer. To explain the issue, let's suppose that you have the following (very small) matrix of counts:

       sample1 sample2
gene1        1       2
gene2       10      20
gene3       15      32
gene4      100     198

As you can see, sample2 has roughly twice as many reads for each gene as sample1, so any normalization is going to cut all of its counts in half (or the reverse for sample1, or perhaps a mixture of the two...depending on the method). That's usually fine, since the assumption with all common normalization methods is that differences like this are simply due to differences in sequencing depth between the samples. But what happens if sample2 actually has twice as much RNA as sample1? Well, then you've just normalized away all of the differences between the samples (I encourage readers to go through the literature on Myc, since that's a very nice example of this, where a mostly-global shift in expression led to apparent down-regulation of genes with constant expression). Such a situation doesn't occur very commonly in bulk RNAseq, but when it does you need to use something like spike-ins (dosed to match cell number or DNA quantity) to maintain the global difference.

If this is only rarely an issue with bulk RNAseq, then why is it more likely to be an issue in scRNAseq? The answer is that bulk RNAseq benefits from the averaging out of many different cell-types, such that ensemble averages/medians don't tend to change much due to treatments...but this isn't going to be the case at the single-cell level, since cells can vary vastly in size (compare a sperm and an oocyte...every RNA should be differentially expressed between them). This difference in base RNA content between different cell types will be washed away by common normalization methods, which can lead to either losing differences or creating differences in genes where there shouldn't be (again, see examples of this in the Myc literature). Whether this really becomes an issue will depend on what your goals are with your data and how much of a global expression difference you actually have and the type of analysis you want to perform.

At the end of the day, just remember that common RNAseq normalization methods assume that (A) there's no global expression change between samples and that (B) if there are very many DE genes that the direction of change is not heavily skewed in one direction. If either of these assumptions are violated then you can run into issues.


Getting sequences of one transcript per gene out of annotation and genome I got a genome (data/genome/genome.fasta) and braker-based genome annotation (data/genome/annotation.gff3), now I would like to get sequences of one transcript per gene (let's say the longest one).

Cufflinks has a gffread program to get the transcriptome out of the genome and the annotation. Something like

gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa

Will do the job. However, it is pulling multiple transcripts per gene. I read the help page twice, found parameters that could potentially do what I would like to (collapsing alternative transcripts into single one), but even when I specify my candidate parameters -M -K I get still .t1 and .t2 sequences for some genes.

The task is surprisingly simple, and I got already a custom script solution, but it just feels wrong. There must be a better way!

What is a reasonable/elegant way to get nucleotide sequences of genes out of a genome and annotation? Am I missing something with cufflinks?

Related, not exactly identical:

How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms (dealing more with terminology of Trinity)
When I used Kallisto/Salmon for differential expression analysis, I collapsed transcripts into genes while importing to R via tximport.

Here is the custom python script solution:

#!/usr/bin/env python3

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt->gene
# jg25690   jg25690.t1
# jg25690   jg25690.t2
# jg25691   jg25691.t1
transcripts2genes = dict()
gene2transcript_mapfilename = 'data/genome/transcripts2genes.map'
with open(gene2transcript_mapfilename, 'r') as g2t_file:
    for line in g2t_file:
        gene, trainscript = line.rstrip('\n').split('\t')
        transcripts2genes[trainscript] = gene

# construct a dictionary gene->list of transcripts
# this fasta file is output of default gffread -w
genes2transcript_seuqneces = defaultdict(list)
transcripts_fasta_filename = 'data/genome/transcripts.fasta'
transcripts_fasta_file = SeqIO.parse(transcripts_fasta_filename, 'fasta')
for seq_record in transcripts_fasta_file:
    gene = transcripts2genes[seq_record.name]
    genes2transcript_seuqneces[gene].append(seq_record)

genes_fasta_filename = 'data/genome/genes.fasta'
with open(genes_fasta_filename, 'w') as out_file:
    # iterate through genes
    for gene in genes2transcript_seuqneces.keys():
        # get all the associated transcripts
        transcripts = genes2transcript_seuqneces[gene]
        # find the longest, by default the one with index 0
        lengest = 0
        for index, tr in enumerate(transcripts):
            if len(transcripts[lengest].seq) > len(tr.seq):
                lengest = index
        # this would be nicer with Bio interface (with wrapping of sequence lines)
        # but I was lazy and blast does not care, only thing I do is converting sequence to uppercase (not even sure if necessary)
        fasta_string = '>' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)


What are the available cloud computing services for bioinformatics? I'm looking for cloud computing services that can be used for doing bioinformatics. An example I found is InsideDNA and there is Amazon of course.
A little description of these would be appreciated.

I have trialed the free version of InsideDNA, and these were my notes:

Cost: \$225/month for a team of 5 with 50TB storage or \$45/month with 10TB storage for individuals (assuming 6 month package: https://insidedna.me/pricing).
Software installed: Around 600 bioinformatic tools available and standard command line tools; some popular tools missing (like CD-HIT), but should be possible to install on request.
Jobs: Maximum of 32 CPUs and 208 RAM per job submission. Test jobs generally worked, although a larger test job failed.
Other points: Command line was sometimes slow, wget queries were slow, and scp was blocked. However, these may be resolvable issues.

Overall, I felt InsideDNA could be useful for groups without their own computational infrastructure and could be used for easily sharing resources between groups. The packages on offer seem not expensive, but I had a few issues, and I don't know how good their sys admin support would be.

I have not used the Amazon service, so can't comment beyond the details on their website. Also there are a few alternative companies, such as Genestack and DNAnexus, but I haven't directly tested them either.


Clarification regarding (deepTools) bamCoverage output Im using bamCoverage on a file and ive set the bin size to 1 but also normalizing using RPKM for defined regions using a bedgraph file.  The output from bamCoverage however is in binned segments each with a (RPKM) score.  Am I correct in interpreting then that
1 - the output bedgraph tells me the (total) RPKM score for that particular bin?
2 - to get the total RPKM score for the region defined in the bedgraph I need to sum all the RPKM scores given for that particular region?

Yes
The total RPKM for the region in the bedgraph file is the value already present in the bedgraph file. If you want to estimate the RPKM for a larger region composed of multiple bins then you can average the per-bin RPKM values. Note that adjacent bins with the same value are merged by deepTools, so you will need to account for this.

Note that RPKM is NOT a robust normalization and that I should remove it from deepTools so people stop using it.


Cannot blast database created with blastdb_aliastool: BLAST Database error: No alias or index file found I want to merge/concatenate BLAST databases as suggested here.

This is my example script that explains what I'm trying to achieve better than words can:

#!/bin/bash

for file in FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn
do
    echo "## Make DB for $file"
 makeblastdb -in $file -dbtype nucl &>/dev/null
    blastn -db $file  -query query_single_nucl.fasta -outfmt 6
    echo
done

echo "## Create testdb.nal..."
blastdb_aliastool -dblist "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -dbtype nucl -title test -out testdb

echo
echo "## Test testdb.nal..."
blastn -db testdb.nal -query query_single_nucl.fasta -outfmt 6

This is the output:

## Make DB for FAM1079.ffn
test_query  FAM1079_00001   100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM3228.ffn
test_query  FAM3228_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM6161.ffn
test_query  FAM6161_000001  100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Make DB for FAM19036.ffn
test_query  FAM19036p_000001    100.000 1266    0   0   1   1266    1   1266    0.0 2338

## Create testdb.nal...
Created nucleotide BLAST (alias) database testdb with 12692 sequences

## Test testdb.nal...
BLAST Database error: No alias or index file found for nucleotide database [testdb.nal] in search path [/path/to/files::]

While the individual databases work, the concatenated databases don't.

The problem was that the database must be specified without suffix, i.e. blastn -db testdb as opposed to blastn -db testdb.nal.

By the way, it produces exactly the same output as

blastn -db "FAM1079.ffn FAM3228.ffn FAM6161.ffn FAM19036.ffn" -query query_single_nucl.fasta -outfmt 6

and the performance is not better.


Albacore wheel not supported on this platform (Mac OS X El Capitan) I download the Albacore wheel from Oxford Nanopore.

I activated my conda environment that has Python 3 as the default python version.

I tried to install Albacore with

pip3 install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

but got the following error:

ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl is not a supported wheel on this platform.
You are using pip version 8.1.2, however version 9.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.`

I then upgraded my pip with:

pip install --upgrade pip
Requirement already up-to-date: pip in /Users/.../anaconda3/envs/python3env/lib/python3.6/site-packages`

So it seems that my pip was already up to date. I do not understand why it complains about my pip if I am using pip3 anyways. I tried to install Albacore again but I got the same message as before.

Can anyone help me install Albacore on my Mac El Capitan laptop, please?

Update based on Wouter's suggestion:

python -m pip install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Processing ./ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl
Requirement already satisfied: python-dateutil in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: setuptools in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Requirement already satisfied: h5py in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting progressbar33 (from ont-albacore==2.1.10)
Requirement already satisfied: numpy>=1.12.1 in /Users/userName/anaconda3/lib/python3.6/site-packages (from ont-albacore==2.1.10)
Collecting ont-fast5-api>=0.4.1 (from ont-albacore==2.1.10)
  Using cached ont_fast5_api-0.4.1-py2.py3-none-any.whl
Requirement already satisfied: six>=1.5 in /Users/userName/anaconda3/lib/python3.6/site-packages (from python-dateutil->ont-albacore==2.1.10)
Installing collected packages: progressbar33, ont-fast5-api, ont-albacore
Successfully installed ont-albacore-2.1.10 ont-fast5-api-0.4.1 progressbar33-2.4

At first sight, it looks like Albacore was successfully installed but if I run: read_fast5_basecaller.py --help

I get:

Traceback (most recent call last):
  File "/Users/userName/anaconda3/bin/read_fast5_basecaller.py", line 69, in <module>
    from albacore.input_utils import list_input_files
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/input_utils.py", line 7, in <module>
    from albacore import log_utils
  File "/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/log_utils.py", line 5, in <module>
    from albacore.pipeline_core import get_debug_level
ImportError: dlopen(/Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so, 2): Library not loaded: /usr/local/opt/python3/Frameworks/Python.framework/Versions/3.6/Python
  Referenced from: /Users/userName/anaconda3/lib/python3.6/site-packages/albacore/pipeline_core.cpython-36m-darwin.so
  Reason: image not found

Is there a library missing?

Not sure about conda as I don't use it, but I can highly recommend pipenv which is the python.org recommended package management system for python.

I have tested the following on my Mac and it worked like a dream

mkdir albacore && cd albacore
mv ../ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl .

# install pipenv if you don't have it
brew install pipenv  # or pip3 install pipenv

# the following will create a virtual environment and install albacore into it
pipenv install ont_albacore-2.1.10-cp36-cp36m-macosx_10_11_x86_64.whl

# activate the virtual environment
pipenv shell

read_fast5_basecaller.py -h

And that should bring up the albacore help menu and you're ready to go.


Rdkit.Chem.Draw.savefig() creates empty figures I can draw a molecule in a jupyter notebook with RDKit. Though, when I use the savefig() method only a white canvas is saved. No molecule is saved.

%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')
fig.savefig('mol.jpeg')

Any idea what goes wrong here?

For some reason you have to set bbox_inches to tight. Without that rdkit does not fit the molecule into the canvas and might plot an empty figure.

fig.savefig('mol.jpeg', bbox_inches='tight')


SnpEff gene annotation with @ such as PCDHB@ I'm using snpEff to annotate some VCFs and getting some spurious annotations that don't have transcript ids. Some of these annotations have a gene name followed by ampersand such as PCDHB@ and I'm not finding any explanation in the docs.

Any idea what this refers to? Entire gene family?

Thanks.

The @ at the end of a name denotes a "gene cluster in chromosome region". This is mentioned in the official HUGO nomenclature guidelines in the "Characters reserved for specific usage" section:

It's also mentioned in the GeneCards documentation:


Identifying Indels from Chromatograms I have around 100 chromatograms (.ab1 files) from Sanger sequencing a genome at loci believed to have an indel.

I'm new to interpreting this kind of data in general, but I've read a bit on the general idea—mostly in guides like this. What I'm not able to figure out from the resources I've checked is how to identify insertions and deletions in a chromatogram.

I'm most concerned with  heterozygous indels, and these seem to be simpler. (And since my samples are (supposed to be) larger indels that are more rare and more likely to be het, I think.) A heterozygous indel, according to my understanding, would like something like consecutive heterozygous SNPs—overlapping peaks, maybe not with the same height, but with the same horizontal placement—until the end of the sample.

In other words, something like the top track in the open window in this image.

(Source: CodonCode.)

(I'm trying to do this manually, and not with packages like this for several reasons, including: I want to see how this is actually done, I want to do it as thoroughly as possible—I will probably check my results with an automated package—and I was asked to do it like this.)

How can evidence of (het) indels like this be distinguished from noise in the chromatogram? My guess is it includes things like

the overlapping peaks last until the end of the sample (and don't start at the beginning, but somewhere in the middle); and
the peaks are still at the same horizontal positions; i.e., one is not shifted right or left relative to the other?

Is this correct?

With my own chromatograms, I have a lot of doubt about distinguishing noise from indel. Here's an example, visualized with 4Peaks. (Full chromatogram here.)

On the left, the chromatogram appears high quality. There are sharp, regularly spaced, single peaks. On the right, there appear to be two sets of peaks lined up with each other. Does this constitute evidence for an indel? Does the homopolymer region—I've noticed in many of these samples, the putative indel seems to start in a homopolymer region—cast doubt on this?

I've also noticed that in some situations like this, there appear to be not two but three peaks.

Is one of these just noise and can be disregarded? Or does this cast doubt on the validity of the data—which is supposed to be coming from the diploid human genome?

If you want to call heterozygous SNPs and INDELs, you may try polyphred. I believe there are a couple of other similar tools, but I have forgotten their names.

  How can evidence of (het) indels like this be distinguished from noise in the chromatogram?

Your two points are correct. Basically, at some position, you start to see two streaks of overlapping peaks. Once identified, the two streaks can be aligned to each other if you shift one of them. This is a strong constraint, but it is hard to tell by eye because it is not obvious how to assign peaks to streaks without running an automatic tool. Doing this is even more difficult given noisy traces like your last figure. I think you'd better run polyphred or similar.

  the putative indel seems to start in a homopolymer region—cast doubt on this?

The INDEL mutation rate in long homopolymer runs is orders of magnitude higher than elsewhere in the genome. This has two implications: 1) you expect to see more homopolymer indels; 2) PCR polymerase is more likely to introduce errors in homopolymers––a higher fraction of homopolymer indels tend to be wrong.


How can I find mutations associated with disease in human histone residues? I would like to look if there are mutations in residues of human histones associated with any disease. For instance, if a mutation in residue K6 (lysine 6) of histone H2A1A is associated with any human disease, (with PUBMED evidence, preferably), but in a systematic way (for all types of histones and all residues).

The purpose of this is to get some candidate residues to study. I have checked the OMIM database but I didn't find anything relevant. It seems that there is no systematic way to connect specific residues to diseases and search for them.

Is there any web database or API available for this purpose?

Ensembl contains this information: When you go to the “phenotype” menu item of a given gene, you will see a list of variants (potentially after clicking on “ALL associated variants”) with their associated genomic position.

You can subtract the transcript’s start position from that position to find out the residue (of the translated sequence, and be careful to subtract the lengths of all preceding introns!). Unfortunately this seems to require some manual work.

Alternatively you can query the individual SNPs, which will tell you their position in all transcripts (in this case: Leu 298/68/244):


Seurat CCA, AlignSubspace, Argument dist.method is only useful with multivariate timeseries I am trying to combine two datasets to be able to compare the cells in them and using Seurat CCA for it:

https://satijalab.org/seurat/Seurat_AlignmentTutorial.html

When I am running AlignSubspace function I am getting a warning for each single dimension that I selected it to be run on:

seurat_union <- AlignSubspace(object = seurat_union, reduction.type =
    "cca", grouping.var = "protocol", dims.align = 1:13)

So, I am getting 13 warnings:

> warnings()
Warning messages:
1: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
2: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
3: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
4: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
5: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
6: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
7: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
8: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
9: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
10: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
11: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
12: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries
13: In dtw(x = align.1, y = align.2, keep.internals = TRUE,  ... :
Argument dist.method is only useful with multivariate timeseries

What is this warning telling me? What should I pay attention to with my datasets?

I don't have enough reputation to comment so I will briefly answer here

This warning is fixed in Seurat 2.3.2. Without upgrading you can ignore this message. It doesn't mean anything for the actual analysis performed so your results won't be flawed.


How to filter the PDB for single-domain proteins with full length 3D-structures solved I am trying to reproduce a machine learning model that has been developed here. As one of the datasets, they use single-domain proteins with full length 3D-structures solved. Since I don't have a background in biology, I am not sure how to filter the PDB for these kind of proteins. Any help would be appreciated.

If you want to reproduce the model in that paper, you should email the corresponding author asking for the exact set of PDB IDs/chain IDs/residue ranges they used. There is too much flexibility in the criteria as given to get the exact set any other way.

If you want to assemble a similar set for your own work, the RCSB PDB advanced search is very powerful and would let you filter e.g. by max sequence length or for only monomers.

However the PDB doesn't deal with domains (regions that evolve and fold independently) directly. A PDB chain may contain part of a domain, one domain, or multiple domains.

If you want a dataset of protein domains you could use CATH, SCOP or ECOD.

As a side issue, please remember to remove overlap of homologous proteins between the training and test sets. People entering the field often neglect to do that. See for example this comment piece.


Running GSEA or comparable analysis with multiple variable interaction? I am working with an RNA-seq dataset generated from a complex experimental design. Our main question is focused on the interaction between two variables, one (INJ) with two conditions (L,S) and the other (SOC) with four conditions (A,B,C,D). We also have a third variable for tissue type (TIS) from 3 tissues (H,N,T), although we are not interested in this variable for this particular question.
We were able extract DE genes underlying the interaction of interest in DESeq2 using an LRT test. We also ran WGCNA to get gene clusters and looked at trait correlations. Using these WGNCA clusters and VST normalized count data, I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables.
What I've tried so far was generating a pseudo-variable that is a combination of INJ and SOC with 8 separate conditions for running in GSEA, but this doesn't quite get at the interaction between variables since we really want to see if the effect of INJ differs between SOC conditions.
One thing I was thinking was running a GSEA-like analysis on log2 ratios between INJ states for each SOC condition, but these would have to be summary statistics for groups of samples and I would lose the variance of expression within groups (maybe that doesn't matter though for exploratory analysis?). This would also reduce my "sample" count to 1 per condition, which isn't ideal for this type of analysis and I would probably have to look at differences between gene sets instead of phenotypes for more robust stats.
Any thoughts on ways to do this? I don't need to use GSEA if there is another analysis I could run on gene sets that allows for more complex experimental models. Thanks in advance.

I'm not sure I understand "I would like to run GSEA or a similar analysis to find WGCNCA clusters that differ based on the interaction between the two main variables". I would run an association analysis (regression) for module eigengenes and the appropriate interaction term. The (most) significant modules are your candidates. This step is simply regression (in the simplest case a linear model should be sufficient), no GSEA needed. You can then subject the identified modules to GSEA or similar to get an idea of their biological function.


How to run same command on multiple files? Hi im doing Variant calling on fastq files for which i have 4000 fq files and the variant calling are done in different 9 steps. Each step generate different files that are output to next steps.
So can anybody help me with this to do it in a loop it pick up names from previous steps for all files. Step 1 command is here:
bwa mem -M -R '@RG\tID:ERR040140\tLB:2517707\tPL:ILLUMINA\tPM:HiSeq2000\tSM:ERR040140' ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta ERR040140_* > ERR040140.sam

So i made an excel file to break the command in columns so that i can change the file names keeping remaining command same just like this:

So changes are occurring 5th, 7th and 9th column. which are:
tID:ERR040140
tLB:2517707
tSM:ERR040140
ERR040140_*
ERR040140.sam

Can anyone help me with loop or something?
Thank you!

Another solution for repeating one command with different parameters is to use parallel.

Create a tab-delimited file samples.txt, which contains the sample names and the information for your read groups, e.g.:

Sample1 Lib1
Sample2 Lib2
Sample3 Lib3

Then you can use parallel like this:

parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:{1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta {1}_* > {1}.sam' :::: samples.txt

For each line in samples.txt a new process is started, replacing {1} and {2} with the values in column 1 and 2 of this line.

fin swimmer


Building RNAz on MacOS Having trouble building RNAz 2.1 on MacOS. Following instructions at https://www.tbi.univie.ac.at/software/RNAz/

tar -xzf RNAz-2.1-tar.gz
cd RNAz-2.1
./configure
make

make fails with

/usr/bin/clang++ -o RNAz RNAz.o rnaz_utils.o svm_helper.o zscore.o cmdline.o strand.o svm.o -L../librna -lRNA -lm
Undefined symbols for architecture x86_64:
  "_HairpinE", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
  "_LoopEnergy", referenced from:
      _alifold in libRNA.a(alifold.o)
      _fill_arrays in libRNA.a(fold.o)
      _backtrack in libRNA.a(fold.o)
      _stack_energy in libRNA.a(fold.o)
ld: symbol(s) not found for architecture x86_64

Thinking that compiling for 32-bit might work I tried setting

CXXFLAGS = -m32
LDFLAGS = -m32

without success (same result). Same result with gmake.

Configure
Configure checks for gcc

...
checking for gcc... gcc
checking whether the C compiler works... yes
...

Makefile

Makefile specifies the following:

CC = gcc
CCDEPMODE = depmode=gcc3
CFLAGS = -g -O2
CPP = gcc -E
CPPFLAGS =
CXX = g++
CXXDEPMODE = depmode=gcc3
CXXFLAGS = -g -O2

Compiler

Apple LLVM version 10.0.0 (clang-1000.11.45.5)

Target: x86_64-apple-darwin18.2.0

Thread model: posix

InstalledDir: /Applications/Xcode.app/Contents/Developer/Toolchains/XcodeDefault.xctoolchain/usr/bin

Has anyone built this successfully on MacOS? Or perhaps suggestions as to how I might build this successfully? I'm at a loss.

Compile with gcc or remove inline from librna/fold.c:776 and librna/fold.c:846.


What assembler is appropriate for High-Fidelity PacBio reads What assembler is appropriate for High-Fidelity PacBio reads?

For example, canu is good for high-error PacBio reads. But what algorithm to use for HiFi reads? Would it be OK to use canu without the error-correction step?

As recently tweeted by A. Phillippy, Canu v1.9 now supports HiFi reads.

Therefore you should just read their manual and look for the implementation

:)


How I can run this code on my files? I am annotating some .txt files by Annovar software by this code

nnovar]$ module load annovar/2016Feb01

[cyan01 annovar]$ table_annovar.pl

But I really got confused to do that by awk for all of my files in the path; Can somebody help me?

Many thanks

  EDITED

I found this but not working

module load annovar/20

for file in *.vcf; do
    table_annovar.pl example/a.txt humandb/ -buildver hg19 -out file -remove -protocol refGene,cytoBand,exac03,avsnp147,dbnsfp30a -operation gx,r,f,f,f "$f" > "${f%.*}.txt"
done;

$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f.

EDIT AS THE OP HAS EDITED THE QUESTION:

First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the right direction.

In you case, the problem is that you are trying to use a variable that you did not declare in the first place. It is explained above in my original answer but here I go again:

With this bit for file in *.vcf;, you are creating a variable file that you will using throughout the iteration. But then within the iteration, the next line (... -operation gx,r,f,f,f "$f" > ..., you try to use another variable file, that has not been declared yet (as far as I can see from your code bit).


10bp sliding window quality score I am having trouble creating a 10 basepair sliding window that goes through seq and trim them all 10 basepairs that have an average quality score (qual) less than 15. I am trying to approach this by reading in the file and seperating the sequences and quality scores, then read in the 10 basepair window in the sequence.
#!/usr/bin/envpython
with open('fastaQ.fq') as f:
    lines = [line.rstrip() for line in f]
    seq= []
    qual= []
    tenFrame= []

    #Seperates the sequences and the quality score of the file
    for i in range(len(lines)):
        if lines[i].startswith('@'):
            seq.append(lines[i+1])
            qual.append(lines[i+3])
            i=i+4

    #
    for i in range(0,len(seq)):
        tenFrame.append(seq[i][i:i+10])
        if len(tenFrame)%10 == 0:
            print(tenFrame)```

I'm not sure what you mean by "trim". The code below will discard any reads that have a 10bp window with a mean phred score below your cutoff. I'll ammend the post to also include a solution that "trims" basepairs on both ends until the window achieves an average score of above the cutoff. It would be ideal if you could provide a minimal reproducible example that shows what your code does and what you'd like it to do.
As far as the code goes, I'd first suggest using a library like Biopython to parse the Fastq file. It will make your life a lot easier. For your code, it looks like seq actually contains all the sequences, so len(seq) is the number of reads in the file. However, when you iterate over i in range(len(seq)), it looks like you're using i to both index which read to iterate over as well as get the 10 base pair window. You need an additional index variable to obtain the window:
for window_start in range(len(seq[i]) - window_len) + 1)
    window = seq[i][window_start:window_start + window_len]

Additionally, some pythonic tips for the future:

If a range starts from 0, don't specify it as an arugment i.e. use range(x) as opposed to range(0, x)
Try to avoid anything along the lines of range(len(x)). In this case, you could just iterate over the sequences directly: for s in seq: or, if you wanted the index, you could enumerate the sequences via for seq_idx, s in enumerate(seq)
Use good variable names. seq implies a single sequence, yet is a list of sequences. Same with qual

Here is some code that will retain only sequences that have all windows with a passing mean phred score
from Bio import SeqIO

read_file = "/file/here.fq"
w = 10
cutoff = 15
valid_seqs = []
with open(fq_file) as fq_fd:
    seqs = SeqIO.parse(fq_fd, "fastq")
    for seq in seqs:
        scores = seq.letter_annotations["phred_quality"]
        mean_scores = (np.mean(scores[i:i+w]) for i in range(len(scores) - w))
        if any(mean_score < cutoff for mean_score in mean_scores):
            print(f"{seq.id} does not pass the filter!")
            continue
        valid_seqs.append(seq)


Getting data from fastq by generator I have a task in a training that I have to read and filter the 'good' reads of big fastq files.
I downsampled, got the code working, saving in a python dictionary. But turns out the original files are huge and I rewrite the code to give a generator. It did work for the down-sampled sample. But I was wondering if it's a good idea to get out all the data and filtering in a dictionary. Does anybody here a better idea? I am asking because I am learning python on my own and I am not sure if it makes sense.

I got some ideas from a code in Biostar:

import sys
import gzip

filename = sys.argv[1]

def parsing_fastq_files(filename):

    with gzip.open(filename, "rb") as infile:
        count_lines = 0
            for line in infile:
            line = line.decode()
            if count_lines % 4 == 0:
                ids = line[1:].strip()
                yield ids
            if count_lines == 1:
                reads = line.rstrip()
                yield reads
        count_lines += 1

total_reads = parsing_fastq_files(filename)
print(next(total_reads))
print(next(total_reads))

I now need to figure out to get the data filtered by using 'if value.endswith('expression'):' but if I use a dict for example, but that's my doubt because of the number of keys and values.

Have a function that returns a single-end read (or pair) and another that does any relevant filtering on that. Wrap both of them in another function that checks if you've hit the end of the file. That's both simple, will handle any size of file, and allows single and paired-end reads.


How to trigger next process on completion of previous process? Simple question but I am stumped. (I am still learning Nextflow)
I have a process which takes as input a channel of files and processes each file separately.
There is a second process which needs to wait until all files are complete, so it can merge the outputs into a single collated file.
How can I do this? I have tried subscribing to "onComplete" event to no avail as I cannot specify a process as the function to trigger.

files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    ????

    output:
    ????

    script:
    """
    echo "How to trigger this process once, after all individual files are processed?"
    """

}

??? Want to do something like this
result.subscribe onComplete: merge_results

```

So from this page...
https://groups.google.com/g/nextflow/c/jt77_-uApMs
It is the "scatter-gather" pattern. I think the answer is to collect() the channel as follows - interested to hear if there is a better way
files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    val go from result.collect()

    script:
    """
    echo "Triggered once after all files complete!"
    """

}


BWA doesn't recognize HiC options I'm trying to run bwa on HiC data. This is the command I'm giving:
bwa mem -5SP -T0 -t4 -R '@RG\tID:\tSM:\tLB:\tPL:ILLUMINA\tPU:none' REFERENCE.fa PAIR1.fq PAIR2.fq

and here's the error message:
mem: illegal option -- 5

Supposedly, the 5SP option is for HiC data (cite), but my BWA isn't recognizing it. Anybody have any ideas? My initial thought is to update BWA to the newest version, but I don't want to mess with our replicatibility since I am new to the lab and don't want to be updating software when I don't need to. The server is a Mac Pro and bwa version is 0.7.12-r1039. Thanks.

Your link with those command line options is 1 year old.  Looks like that version of bwa is 6 years old.  You probably should update.


GATK CombineVariants complains the contig order in the VCF files I have called variants on two strains of C. elegans separately. I now want to merge the VCF files into one using the following code:

Create a sequence dictionary of the reference sequence
Sort the VCF files with Picard
Merge the sorted VCF files using GATK

picard CreateSequenceDictionary \
    REFERENCE=c_elegans.PRJNA13758.WS263.genomic.fa \
    OUTPUT=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain1.vcf \
    OUTPUT=strain1sorted.vcf \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

picard SortVcf INPUT=strain2.vcf \
    OUTPUT=strain2sorted.vcf  \
    SEQUENCE_DICTIONARY=c_elegans.PRJNA13758.WS263.genomic.dict

GATK --analysis_type CombineVariants \
    -R c_elegans.PRJNA13758.WS263.genomic.fa \
    --variant strain1sorted.vcf \
    --variant strain2sorted.vcf \
    -o all.vcf \
    -genotypeMergeOptions UNIQUIFY

The last command gives me the following error message:

ERROR MESSAGE: Input files variant and reference have incompatible contigs. Please see https://www.broadinstitute.org/gatk/guide/article?id=63 for more information. Error details: The contig order in variant and reference is not the same; to fix this please see: (https://www.broadinstitute.org/gatk/guide/article?id=1328),  which describes reordering contigs in BAM and VCF files..
##### ERROR   variant contigs = [I, II, III, IV, MtDNA, V, X]
##### ERROR   reference contigs = [I, II, III, IV, V, X, MtDNA]

But I have sorted the VCF files using Picard, so I don't know what else to do.

Your help is appreciated.

Update: I asked this question on the GATK forum on 4th May 2018 (today).

Answered in the GATK forum

  Try regenerating the index files for your VCFs. Picard SortVcf doesn't do it for you iirc, so when GATK looks at the index files (before opening the VCFs themselves) it sees the old order and complains.


tools to reconcile experimental transcripts with reference annotation Looking for tools to reconcile alignment file of experimental transcripts mapped to genome (SAM/BAM) with the reference transcriptome annotation (GTF) from Ensembl (organism: D. melanogaster).

The aim would be to check which transcripts reported in the alignment file are also reported in the reference annotation, generating a new annotation (e.g., GTF) including the new transcripts as seen in the alignment file (experimental transcripts to genome).

Any idea?

-- Edit:

Up to now, I tried an in-house script that does the following:

generate a BED file from the BAM file with the exon coordinates (start/end) per transcript - using BEDOPS v2p4p27
convert the reference GTF to BED to retrieve the coordinates per annotated feature
compare the experimental vs the reference BED intervals
infer the suite of annotated exons per experimental transcript, thus, ideally, the isoform it corresponds to

However, with this simple approach, I only get around 20% of the experimental transcripts that can be identified in all their length using reference annotation. This might be due either to a too simple approach (e.g. if the start/end of the exon in the experimental transcript is shifted by 1 base as compared to the reference, it won't be detected), or to a real biological signal (since I am also using data from long-read sequencing technologies).

I've never tried this myself, so I don't know how easy this is...

One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You should then be able to get the information you want from the transfrag codes.


You are almost there except for a few tricks:

If you want columns to be ordered using the values in Cell 4, you will need to provide a custom column order either with the column_order argument of the Heatmap() function or by providing ordered data. Below is the code to perform the latter. Most of the code is yours, I just tweaked it a little.

library(tibble)
library(dplyr)
library(ComplexHeatmap)
library(RColorBrewer)
library(circlize)

# Since you want to use cell4 as reference
cell4 <- as.data.frame(t(data["Cell4",]))

df <- cbind(cell4,
            Group = as.character(ann$Group))

# ordering by Group and by expression value, minus denotes decreasing order
df <- df[with(df, order(Group,-Cell4)),]

# to make it compatible to your code in your post
df <- rownames_to_column(df, "Samples")

#Set annotation
ann <- data.frame(df$Group)
colnames(ann) <- c("Group")
colours <- list("Group" = c("A" = "orange", "B" = "darkgreen"))
colAnn <-
  HeatmapAnnotation(
    df = ann,
    which = "col",
    col = colours,
    annotation_width = unit(c(1, 4), "cm"),
    gap = unit(1, "mm")
  )

# rownames of the annotation df and the colnames of the data should be in
# the same order
data <- data[, df$Samples]

myCol <-
  colorRampPalette(c("#ffffff", "#0030dd", "#3259e3", "#7f97ee"))(100)
myBreaks <- seq(0, 1, length = 100)
hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = TRUE, # you asked your cells to be ordered
    cluster_columns = FALSE, # column order is the same as the provided data
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right"

)

The resulting heatmap, still not quite what you want:

You will realize that the values for Cell 4 are ordered from high to low for each of the two groups. To make Cell 4 be drawn at the top, you will need to reorder the rows. To do so, you will need to extract the row orders with the row_order() function:

> row_order(hmap)
 [1]  3  2 14  4 12 11  9 10  7  1  5 13  6  8

Now you can re-draw the heatmap after altering the row order somehow. For this sample case, doing it manually.

# promoting Cell 4 to the first position
modified_row_order <- c(4, 3, 2, 14, 12, 11,  9, 10,  7,  1,  5, 13,  6,  8)

hmap <-
  Heatmap(
    as.matrix(data),
    name = "xCell scores",
    col = colorRamp2(myBreaks, myCol),
    width = unit(100, "mm"),
    show_row_names = TRUE,
    show_column_names = FALSE,
    cluster_rows = FALSE, # set to FALSE as a custom row order is supplied
    cluster_columns = FALSE,
    show_column_dend = FALSE,
    show_row_dend = FALSE,
    row_dend_reorder = FALSE,
    column_dend_reorder = FALSE,
    clustering_method_rows = "ward.D2",
    top_annotation = colAnn,
    column_order = NULL,
    row_order = modified_row_order # argument to supply custom row order
  )
draw(hmap,
     heatmap_legend_side = "right",
     annotation_legend_side = "right")

If I have misunderstood your question and you do not want columns to be ordered according to expression values in Cell4, forget about the first part with the ordered annotation data frame.


Where to Download Cancer Raw Reads (fastq)? Does anyone know where I can download cancer raw reads (fastq files), tumor and Germline for non-humans?

I wanted to make a study with human data but I don't have the control access to download raw reads.

You can try to browse SRA (sequence read archive) at https://www.ncbi.nlm.nih.gov/sra.
To download fastq files from SRA, you'll need their sratoolkit: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software.

The easy command to download fastq files by their accession number is :

fastq-dump --split-files accession_number


Selecting 65000 SNPs where AF is close to 0.5 in all or most populations I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.

My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.

I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that:
(1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible.
(2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.

My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:

gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /\;AF=0\.[34567]/) { if ($_ =~ /\w+\t\d+\t\S+\t\w\t\w\t/) {print $_} }'

Any recommendations?

Just use bcftools view for filtering:

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf

To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.

$ bcftools view -h input.vcf.gz|wc -l
255

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf


What origins for virus HIV1? What origins for virus HIV1?
This is the alignment of a small part of the POL protein of the HIV1 virus GenBank KU749412.1, with the POL protein of the Visna virus GenBank L06906.1
IKKKDSTKWRKLVDFRELNKRTQDFWEVQLGIPHPAGLKKKKSVTVLDVGDAYFSVPLDESFRKYTAFTI--PSTNNETPGIRYQYNVLPQGWKGSPAIFQSSMTKILEPFRAQNPEIVIYQYMDDLYVGSDLEIGQHRAKIEELREHLLRWGFTTPDKKHQKEPPF--------LWMGYELHPDKWTVQPIQLPE--KDSWTVNDIQKLV
|||| | ||| | ||||||| | |  | ||| ||| ||  || || || |||||  || |  | || |    |  ||  |  ||   ||||||| |||  |  | |||       | |    |||| | ||||  |      || |                                | |||| ||  |   |||      | |  ||||
IKKK-SGKWRMLIDFRELNKQTEDLAEAQLGLPHPGGLQRKKHVTILDIGDAYFTIPLYEPYRQYTCFTMLSP--NNLGPCVRYYWKVLPQGWKLSPAVYQFTMQKILRGWIEEHPMIQFGIYMDDIYIGSDL--GL-----EEHRGIVNELASYIAQYGFMLPEDKRQEGYPAKWL-GFELHPEKWKFQKHTLPEITEGPITLNKLQKLV

Are these matches,

By chance or,
HIV1 is a direct descendant of the Visna virus or,
they have a shared progenitor?

So it could also be that HIV1 really comes from some kind of experiment involving Visna virus, as I read on the Internet.
The disease caused by the Visna virus is (100% deadly see https://pubmed.ncbi.nlm.nih.gov/16425963/ ). Could there be parallels to sars2?
But look at it: the Visna virus is transmitted by air the first symptoms are like pneumonia (or a brain infection) and otherwise remains slowly progressing... like sars2, see long covid...
Symptoms are also felt after many days that the virus has come into action, this seems to me not to happen for coronaviruses or for flu viruses but for HIV and Visna virus
Link related:

https://pubmed.ncbi.nlm.nih.gov/8289690/

https://pubmed.ncbi.nlm.nih.gov/16425963/

(preprint sars2 as retrovirus) https://www.biorxiv.org/content/10.1101/2020.12.12.422516v1

https://www.gao.gov/assets/100/91319.html

They are both lentiviruses and share a distant common ancestor.
HIV-1 and HIV-2 are descendents from simian immunodeficiency virus (SIV). The following tree shows the relationships very clearly, from Wertheim and Worobey (2009) Dating the Age of the SIV Lineages That Gave Rise to HIV-1 and HIV-2 PLoS Comp Biol. here. The evolution of HIV-1 is heavily mixed with chimanpzee (cpz) SIV, and there are regions of highly anomalous homology juxtaposed with regions of notably lower homology.


Dataset showing evolution of gene As a way to exemplify how the hierarchical clustering of particular DNA sequences can be used to derive phylogenetic trees, I'm looking for a (even very small) dataset of the same genes (or other biological sequences) from different organisms (for example animals where the evolutionary tree is known). Do you have any pointers?
Thanks in advance

If you're building a tree of different species, it's probably better to use protein sequences. A common resource is Uniprot.
If you still want DNA sequences you can turn to NCBI Nucleotide.
Some examples that give decent trees: p53, TOP1, HBA.
A common way to construct a phylogenetic tree is to start with a multiple sequence alignment (MSA).


Nanopore data clarification regarding merging of samples We had same set of sample ran on nanopore platform twice.
The issue on hand is.
First run was done it was nearly done but the power went off.
Second run was done completely with enough sequencing depth.
Now the question I ran the analysis with the second set which was complete that set of result I have. But so far I have't ran the first part.
Question1 : if I run the first part do i get major differences? in the result
Question2 : would it be advisable to merge both run since they are now can be called as twi different batches plus would be bring drastic changes in the output ? And in terms of tool for merging I googled which shows poretools.

It is entirely unclear which type of analysis you aim to do. You probably should be more specific when asking questions.
That said, given that the run was going perfectly fine until a power issue stopped it, all data produced should be fine and you can indeed combine those runs.
You probably shouldn't bother with poretools. If you are using fastq files to start your analysis you can just cat the basecalled files together. Also after alignment, you can merge the bam files. If your analysis requires the fast5 files then you can probably just put these in the same folder, but you should tell us which analysis you are doing.


How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms I am working with RNA-seq data without a reference genome/transcriptome and am instead using a Trinity de novo transcriptome assembly.
I analyzed both isoform and gene expression abundances using RSEM.
And got differentially expressed GENES via DESeq2.
In the next step I want to get fasta sequences of differentially expressed genes to use them in the subsequent BlastX/Blast2Go analyses.

Considering the fact that my Trinity created assembly is based on isoform id:

How can I extract the gene sequence out of it? Is that a right thing to combine all the isoforms of a single gene in order to generate the longest possible isoform assembly?

If so, how can I do it to have a union-exon (combined form of all the exons present in the isoforms of an individual gene) sequence used in my Blast analysis?

If not, how should I proceed with the data type available?
Trinity assembled file includes the isoform column as:

TRINITY_DN1780_c0_g1_i2
TRINITY_DN1780_c0_g1_i1
TRINITY_DN1780_c0_g1_i3
TRINITY_DN1780_c0_g1_i4
TRINITY_DN1780_c0_g1_i5

and the DESeq2 produces the output (gene differentials) as:

TRINITY_DN1780_c1_g1

From your question it isn't clear if you are trying to identify the isoforms that have a very high or very low differential expression value, or if you are trying to correctly find an orthlog for your top hits using BLAST.

If your question is about isoform-level quantification, identifying DE at the isoform level is challenging and DESeq2 is only designed for gene-level analyses. That is why it is only outputting genes rather than a specific isoform.


GRCh38 vcf file with common cancer mutations Is there a vcf file on the GRCh38 assembly with common cancer mutations I can download somewhere? Maybe from one of the big international cancer genomics consortia?

By common, I mean whichever mutations have been found recurrent in different types of cancer.

Your question doesn't give enough information for a specific answer but this should do for a start.

Get the VCF file describing all variants in Clinvar from NCBI:

wget ftp://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz

Extract any variants whose VCF line contains the word "cancer":

zgrep -iE '^#|CLNDBN=[^;]*cancer' clinvar.vcf.gz > cancer.vcf

Now, you should probably filter that further based on variant frequency, type of cancer, etc etc but that should be a good start.


Looper: can't specify the individual path for one pipeline run I have multiple samples that I need to run one pipeline on. I am using Looper: https://github.com/pepkit/looper

I have a real trouble specifying a separate location for each pipeline output files. The first time when I did not care about it thinking that Looper would just automatically do that, it lead to multiple files being overwritten. Then, trying to avoid it I tried (and it did not work):

  Adding {sample_name} to the output_dir of the metadata section of project_config.yaml file:

 metadata:
  sample_annotation: sandro_annotation.csv
  output_dir: /scratch/nv4e/sandro_looper_results/{sample_name}
  pipelines_dir: [../pipeline/pipeline_interface.yaml]

  Creating a derived_column in the .csv annotation file:

sample_name,library,sample_name_1_read,sample_name_2_read,output_folder
abeta_24h_1,RNA-seq,read1,read2,out
abeta_24h_2,RNA-seq,read1,read2,out
abeta_24h_3,RNA-seq,read1,read2,out

pipeline_interface.yaml:

...
pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-S": sample_name
      "--input": sample_name_1_read
      "-O": output_folder
...

and project.config.yaml:

...
derived_columns: [sample_name_1_read,sample_name_2_read,output_folder]
data_sources:
  read1: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r1.fq.gz
  read2: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r2.fq.gz
  out: /scratch/nv4e/sandro_looper_results/results_pipeline/{sample_name}
...

In this particular case the parameter -O was formed correctly which I was able to see in sandro_looper_results/submissions/ .sub file, however, it was supplied the second time, apparently from the output_dir section of the project_config.yaml and was overriding the first -O command. So, obviously I tried to delete output_dir line but it is a required setting for Looper, so it can not run without it.

  I have also tried supplying -O command in project_config.yaml in the pipeline_args: section with no success.

I really do not know what I could do at this point. Any suggestions would be greatly appreciated,

Looper does not generate any of these per-sample outputs.

Looper does not have any control over where your pipeline puts its outputs.

The problem is in your pipeline, not in looper. Assuming you're using a pypiper pipeline, I suggest reading the docs; for example, the "your first pipeline" tutorial: http://pypiper.readthedocs.io/en/latest/tutorials-first.html

For real-life example, see this pipeline

outfolder = os.path.abspath(os.path.join(args.output_parent, args.sample_name))
pm = pypiper.PipelineManager(name="RRBS", outfolder=outfolder, args=args, version=__version__)


What is a standard approach to binarize microarray gene expression data? I am using gene expression data from a microarray dataset for the purpose of learning the network structure of genes in the dataset.

Now to do this I need to binarize my data, i.e convert the data for each gene to 0s and 1s. I am new to bioinformatics and am not sure what approach is considered standard for microarray data.

I appreciate your inputs regarding this and I am implementing this in R. So if you know any package that does this please let me know.

Based on the info you provide, ArrayBin R package provides you the necessary tools:

binarize.array() from ArrayBin, allowing:

  Implementation of an adaptive method for binarizing gene expression data on a per-probe basis and demonstrate the superior effectiveness of our method when compared with other, commonly used approaches. This adaptive binarization method can be applied to DNA methylation microarray data,
  which has implications for cross-platform integration, and can reduce batch effects in the data.

In complement, I'd advise you go through some literature and choose the approach which fits you best. Some suggestions:

Computational analysis of microarray data
Binary matrix factorization for analyzing gene expression data
and probably most relevant for you given your application: Multiscale Binarization of Gene Expression Data for Reconstructing Boolean Networks

Maybe you also want to look into more general concepts of transformation prior to microarray data analysis:

Microarray data normalization and transformation
A Tutorial Review of Microarray Data Analysis


Is there an R package that computes homoplasy excess ratios (HER)? I want to calculate the homoplasy excess ratio for a phylogenetic tree, for example homoplasy resulting from recombination.
I am wondering if there are any R packages out there that will do this?
I have looked at phyloseq, Rphylopars, phylocurve, PhylogeneticEM, PhyloinformR, but these packages do not have functions that calculate HER. I searched "homoplasy excess ratio" on bioconductor, but it returned ca. 3000 packages.

Maynard's Smith's homoplasy is what you want and it is available in the recombination detection package (RDP) described here. Its really a Windows package rather than R.


I have 23andme text files and would like to convert to SAM/BAM format I would like to convert 23andme text file to NextGen to BAM file for Yfull.com to read.  It is difficult to get answers on how to convert to SAM/BAM file for the 23andme text file to be converted to SAM/BAM tools

Having done 23andme myself I can tell you that your variant file, which contains SNP genotypes, cannot be converted to a bam file. It does not contain the same information as a bam file. It may be helpful to familiarize yourself with those filetypes and the technologies used to obtain them. A SAM/BAM file contains alignments of reads obtained by sequencing. Your standard* 23andme analysis doesn't have reads, as it uses a SNP array for determining the alleles at predefined positions. Markers which are not included in your 23andme file cannot be added by other means because those were not genotyped by the array. Therefore also Y-DNA cannot use your 23andme file or add additional markers.

*For the sake of completeness: I am aware that under rare circumstances 23andme can perform WGS for additional variant typing, but I do not assume this is the case for OP.


umap and Louvain clustering on normalized data I know that scaled data must be used for PCA for example as it is based on variance maximization.

However I'm wondering if it's the same case for UMAP ? If the data are single-cell RNA seq, after normalization can we do UMAP and Louvain clustering ?

I tried both and get similar results, however the Louvain clustering seems to be more adequate on normalized data than on scaled data.

By adequate I mean the clusters are the same but some are split into two, which makes sens looking at other results (like transcription factor analysis).

If the data are single-cell RNA seq, after normalization can we do
  UMAP and Louvain clustering ?

Typically people run PCA, UMAP and Louvain clustering on the normalised and log-transformed expression counts (but do marker gene and differential expression analysis on the non-normalised values). However, this remains controversial. I recommend reading "Current best practices in single‐cell RNA‐seq analysis: a tutorial"- it's a bit outdated, but still a good reference.

  By adequate I mean the clusters are the same but some are split into
  two, which makes sens looking at other results (like transcription
  factor analysis).

The number of (or size of) clusters should not be used to interpret how good the clustering methodology is, since it is generally dependant on the input parameters, e.g. a higher resolution parameter leads to more clusters. Instead, display the clusters on the UMAP and see if they are continuous 'visually' and check what genes are differentially expressed between clusters and see if they make sense 'biologically'.


How to cluster the human genes by pathways/system-biology/metabolic properties? I would like to make a chord plot from my data. I have a list of genes that according to my experiments are divided into 64 clusters of enrichment pattern.

I would combine my 64 clusters with a second clustering. A clustering of all the human genes into subgroups of "interacting" genes (e.g. both are pathway, co-expressed) will allow me to combine the two clustering data on a chord plot.

In the analogy of the following graph the continents will be my 64 clusters of the genes and the edges/chords will be assigned based on the biological function  clustering. I guess that such a functional clustering can be done ore has been done using platforms such as: KEGG, genescf, or reactome.org.

Is there a biological pathway clustering of genes for lymphoblastoid or breast cell lines?
Are there any suggestions how to build the clustering between genes whose functions are strongly correlated in the cell?

Edit: the question is about getting the clustering from a system-biological perspective not how to plot chord diagrams.

Is there such a metric to quantify how likely two genes are correlated in their enrichment, function etc. For example, using STRING we can see that PIK3CA and PTEN are more co-functioning than PIK3CA and SF3B1. How can I do something like String on genome wide scale.

To measure if two genes are functional similar I developed the BioCor package. It calculates a similarity score between genes by the measuring the amount of pathways  shared. However, it doesn't take into account the strength of the connection between two proteins. You could multiply the similarity by a value between 0 and 1 according to the number of evidences of the relation between the two proteins.

It is optimized for large amounts of data, you could even perform them in parallel, but I doubt you find information about pathways for most of the genes (In my informal test, I think that I could calculate similarities for about 1600 genes not the ~20000 of the human genome).

When I am interested in co-expression data I calculate the co-expression by correlation and add the similarity of the genes calculated with BioCor. The resulting distance matrix can be used for clustering. Maybe the section about using weighted co-expression network analysis and BioCor can be helpful.

To have a "biological pathway clustering of genes for lymphoblastoid or breast cell lines", you would need to have a specific signature/pathway of the cell lines. Currently I am not aware of accepted standard signature database. The single cell atlas could be useful to find expression of your cells, but I don't know if they provide a signature for cell line


NO_COOR reads not in a single block at the end 0 -1 I am doing an RNA seq analysis with mouse samples. I had my reads in several lanes, so I concatenated those fastq reads. Then I've mapped my reads with Hisat2 against the new version of mouse  "GCA_000001635.9_GRCm39_genomic.fna".
Apparently everything works well up to here. Now I'm using samtools 1.9 to convert sam to bam, to sort those bam and create the index, just like this:
samtools view -bS ${file} > ${file}.bam

samtools sort -n ${file} -o ${file}_sorted.bam

samtools index  ${file}_sorted.bam

But when indexing I get this error:

[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 0 -1
samtools index: failed to create index for
"CD4_1A_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...
[E::hts_idx_push] NO_COOR reads not in a single block at the end 23 -1
samtools index: failed to create index for
"CD4_1B_concat.bam_sorted.bam": No such file or directory
[bam_sort_core] merging from 11 files and 1 in-memory blocks...

And more...
Any help would be appreciated!

This is essentially a duplicate of Samtools Index: Chromosome Blocks not Continuous but with different error messages, and a duplicate of Error given while trying to index a BAM file with Samtools Index - NO COOR?.
samtools sort -n ${file} -o ${file}_sorted.bam

These error messages also indicate that the reads are not sorted by coordinate — in particular, that the unplaced unmapped reads with RNAME * are not all together at the end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


ethnicity check either from bam or vcf files What tool could I use to check the ethnicity of a human bam or vcf file? I would like to use the results as a QC check to know whether a given sample or set of samples match the ethnicity information annotated in the meta data or not.

Seen so far:

Similarity:

The examples work, but starting from a GRCh38 .vcf file (from GATK UG) gives me an error:

Read    mytest.vcf   :1 Individual(s)Exception in thread "main" java.lang.NumberFormatException: For input string: ""
        at java.lang.NumberFormatException.forInputString(NumberFormatException.java:65)
        at java.lang.Integer.parseInt(Integer.java:592)
        at java.lang.Integer.valueOf(Integer.java:766)
        at Individual.get_Individual(Individual.java:73)
        at similarity.main(similarity.java:80)

NGSrelate:

It won't work with GATK UG .vcf files straightaway, it needs to reprocess .bam files using the ANGSD pipeline:

VerifyBAMID:

This tool does not give ethnicity checks but only genotype checks. Running all-vs-all vcfs vs bams can create a semi-matrix of similarities.

The main difficulty here is the use of GRCh38.  Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference.  It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the best.

Let's suppose you go with overlapping rsIDs, and a list of those rsIDs, one per line is in 'rsids.txt'.  (Due to strand flips between reference builds, you may want to restrict this list so that all allele codes match, and it's also reasonable to throw out A/T and C/G SNPs.)  Then, the following process would work:

Download ADMIXTURE (https://www.genetics.ucla.edu/software/admixture/download.html ), plink (https://www.cog-genomics.org/plink/1.9/ ), and/or chr1-chr22 of 1000 Genomes phase 3 (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ ) if you don't already have them.
Use plink to extract the overlapping higher-MAF rsIDs from each 1000 Genomes phase 3 VCF ("plink --vcf ... --extract rsids.txt --maf 0.05 --make-bed --out ..."), and then merge the resulting per-chromosome filesets ("plink --merge-list ... --out merged_phase3_subset").
ADMIXTURE prefers a dataset with around 100k variants.  Unless you have a small list of overlapping rsIDs (in which case you had better use another method), you should have more than that.  plink's LD-pruning function is a good way to select a subset for ADMIXTURE's use:
"plink --bfile merged_phase3_subset --indep-pairwise 500kb 1 0.2; plink --bfile merged_phase3_subset --extract plink.prune.in --make-bed --out admixture_data"

Adjust the 0.2 threshold as necessary to keep the right number of variants.

Run ADMIXTURE in unsupervised mode ("admixture admixture_data.bed 5 -j8"; adjust the -j parameter depending on the number of processor cores).  This generates an admixture_data.5.P file with population allele frequencies, and an admixture_data.5.Q file with sample population assignments.  Verify that the .Q file corresponds to the 1000 Genomes phase 3 'superpopulations'.
Convert your data to plink-format if necessary, keeping just the overlapping rsIDs and 'downgrading' to GRCh37 coordinates.  Make sure your variants are sorted in GRCh37-coordinate order, and the allele order also matches that in your reference dataset ("--a2-allele admixture_data.bim 6 2" during your final --make-bed operation will do the trick); otherwise ADMIXTURE won't do the right thing.  Then run ADMIXTURE in projection mode ("cp admixture_data.5.P my_data.5.P.in; admixture -P my_data.bed 5").  my_data.5.Q will then have the ethnicity estimates you're looking for.


How to run MaxQuant in command line mode? MaxQuant is a software package for mass spectroscopy and proteomics. There is a windows version and a linux version. To run on linux you have to use a program that is called mono. I think, it is developed by Microsoft, which I find quite nice from them. I installed mono on Linux and started MaxQuant.

mono MaxQuant.exe

The documentation claims it is for running MaxQuant in commandline mode, though a GUI is started anyway. How can I start the program in commandline mode to perform batch analysis? I think, there is not enough information given on the web-page.

Give

mono MaxQuantCmd.exe [OPTIONS]

a try.

University of Oslo description of how to use the tool in their infrastructure.


Where can I get ensembl Medaka genome for RNAseq I am trying to map a RNA-seq dataset using ensembl genome for medaka fish.

From here http://uswest.ensembl.org/Oryzias_latipes/Info/Index

when I click on a Download DNA sequence (FASTA) I am forwarded here ftp://ftp.ensembl.org/pub/release-98/fasta/oryzias_latipes/dna/

where I have a huge number of .fa.gz files and I am just looking for a single file, that I can use.

While when I am downloading medaka genome from NCBI (and using this technique) http://www.metagenomics.wiki/tools/fastq/ncbi-ftp-genome-download I have my ncbi genome in one file.

How can I get an ensembl genome in one file?

Download the ones without a chromosome number in the filename. Those files are the whole genome.


AutoDock 4 -soft docking How could I turn on/off the smooth potential in AutoDock4? I need to do two docking. The first one is soft docking with smooth potential and the second is re-docking without the smooth potential.

The smooth parameter is an optional keyword and command that changes the smoothing potential. If not provided, the default is 0.5 angstroms by default. Try passing 0 with this keyword in the AutoGrid Parameter File to turn off smoothing potential.

See the AutoDock 4.2.5 documentation for more details.


Random networks in PPi I am working with an in silico model to measure the distances between drugs targets and disease proteins of a disease. I would want to compare these distances with the distances of random networks but I don't know which is the best way to do the comparison. For instance, if I have 300 drug targets and 500 disease proteins mapped in the interactome and I want to measure the distance between them and then compare to random networks with the same features.

How do I generate these random networks? Choosing the same number of drug targets (300) and disease genes (500) randomly in my referral network? Choosing the same number of nodes (300 and 500) and the same number of links?

I'd appreciate your help!

If your goal is only to account for degree distribution (i.e. Matteo Ferla's comment), a simulation approach you could use is:

measure the degree (number of links) of each protein in your PPI
find a large set of proteins that are equal or similar in degree to each of your proteins of interest (you may have to use a range of degrees to get a large enough set for comparison, ideally >=1000). We'll call these "comparable proteins" for each of your proteins of interest.
For each protein of interest, make all of your measurements also for its set of comparable proteins
Compare your measurement (shortest distance?) for the protein of interest to that for its comparable proteins. You could compute Z-scores from the parameters of the comparable protein distribution, or just empirical Monte Carlo p-values from where it falls in the distribution.

This will give you at least a heuristic null distribution for each protein. It should be relatively easy to code up in whatever your favorite language is.

Hope that helps-


How many reads do I need to cover the entire genome? Suppose my genome is 3 million bases and that my reads are 100 nucleotide long. I need to know how many reads I need to cover the entire genome.

I start from using the equation $C = \frac{N \cdot L}{G}$ where C is the coverage, N the number of reads, L the length of a read and G the length of the haploid genome. I also know that the average coverage follows a Poisson distribution, so I need this for C.

On my slides, in fact, I have the following:

$$
C = \frac{N \cdot L}{G} \approx ln \left( \frac{G}{L \cdot \epsilon}\right),
$$

if we assume that this is the average coverage needed to cover the entire genome with probability $1-\epsilon$.

The problem is that I do not understand how to arrive at the expression on the right. I thought that if we want to be sure to have a coverage of at least 1 on average, we would have had to compute $1 - P(X=0) = 1 - e^{-\frac{N \cdot L}{G}} = 0.99$, if we take $\epsilon=0.01$. But I do not get the same result as with the equation above.

Let's see what your slides are claiming $\epsilon$ to be, using:

$$ C = \frac{NL}{G} \approx \ln \left( \frac{G}{L\epsilon} \right) $$

We can rearrange to get

$$ \epsilon \approx \frac{G}{L} \times e^{-C} $$

Now, what does it mean for an entire genome to be covered? Having a clear understanding of this definition is crucial, and your misunderstanding seems to partially stem from this.

I claim it's natural to consider a genome as consisting of $G/L$ bins. In this setting, we have $N$ balls that each randomly goes into a bin. Then, the entire genome is covered when all bins contain at least $1$ ball/read.

From the Poisson distribution, we know that P(a bin has 0 balls) = $e^{-C}$.

Then, we define $\epsilon$ to be the probability that at least one of the $G/L$ bins has 0 balls. Given the way we've modeled the problem, this is exactly equal to

$$ \epsilon = 1 - (1 - e^{-C})^{G/L}. $$

From above, your slides are claiming $ \epsilon \approx \frac{G}{L} \times e^{-C}$. This might be strange, as it appears to be equal to the expected number of bins with 0 balls, instead of the probability that there is at least 1 bin with 0 balls. Furthermore, it's strange because it contains only one exponential, compared to our exact derivation which has a double exponential.

The mismatch in the number of exponential is a hint that another mathematical approximation is at play. We can use a first-order approximation of the binomial series to obtain

$$ (1-e^{-C})^{G/L} \approx 1 - (G/L)e^{-C} $$ which holds particularly when $|e^{-C}| \times \frac{G}{L} << 1$.

In your case, $\frac{G}{L} = 30,000$, so the approximation only really starts getting accurate when $C > 15$ or so.

Following through with our approximation, we get

$$ \epsilon \approx 1 - (1 - (G/L) e^{-C}) $$

$$ \epsilon \approx \frac{G}{L} \times e^{-C}. $$


How to identify genes from a genome assembly of C. Elegans? I have two full genome assemblies for C. Elegans samples collected from two different geographical areas that I found on WormBase. These are in fasta format. I want to go gene-by-gene and compare the nucleotide sequences corresponding to the gene between the two samples. My goal with this is to count how many single-base differences there are in each gene. However, I'm not sure how to match the sequences in the assembly to known reference sequences (by the way is there a good resource to download all of these at once?) for C. Elegans genes.
My thought was to run BLAST between the reference sequence and each genome and find the match with the best score, if any. Is this a reasonable approach, or is there some better way to do this?

You don't want BLAST here, or at least not regular BLASTn. A better method would be to use tBLASTn to map the known C. Elegans proteins to the translated genome. However, I would recommend using a slightly more sophisticated approach:

Collect fasta sequences of all known C. Elegans proteins (amino acid sequences, not nucleotide). This should be easy enough to do on WormBase. For example, the protein sequences for C. elegans strain VC2010, assembly PRJEB28388, can be found here: ftp://ftp.wormbase.org/pub/wormbase/releases/current-production-release/species/c_elegans/PRJEB28388/c_elegans.PRJEB28388.WS279.protein.fa.gz

Use a tool that can model introns and exons to map a protein into a genome. I haven't worked in this field in more than 10 years now, so it is very likely that there are other programs around today, but back in the day, I would do this using either exonerate or genewise. With exonerate, the command would be:
exonerate -m protein2genome -n 1 -t assembly.fasta -q allproteins.pep > out

That will give you a list of best-hits for each of your input proteins and should accurately map the protein to the genome. Once you have that, you can extract the matches and start comparing.


why in RNA seq don't we only use reference transcriptome? I would like to ask why in RNA seq analysis (alignment step) we use sometimes reference genome instead of reference transcriptome? thank you!

Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.
But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.


WGCNA: Problem with selecting soft threshold I have 25 tumor samples with counts data. Initially, I filtered out low expressed genes and then converted counts to varianceStabilizingTransformation using deseq2 package. With this data I started using WGCNA for co-expression network analysis.

For selecting the soft threshold I see very strange plot. R2 cutoff is 0.8 and I see that none of the scale free topology model fit is above that.

Here is the code I used:

df is a dataframe with genes as rows and 25 samples as columns with counts data.

library("DESeq2")
filtered.counts <- df[rowSums(df==0)<5, ]

U3 <- as.matrix(filtered.counts)
vsd <- vst(U3, blind=FALSE)

oed <- vsd

gene.names=rownames(oed)
trans.oed=t(oed)
dim(trans.oed)

n=16462;
datExpr=trans.oed[,1:n]
dim(datExpr)

SubGeneNames=gene.names[1:n]

library(WGCNA)
options(stringsAsFactors = FALSE);
allowWGCNAThreads()

powers = c(c(1:10), seq(from = 12, to=20, by=2));
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
                      powerVector = powers,corFnc = cor,
                      corOptions = list(use = 'p'),networkType = "unsigned")

   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1 0.000273  0.041          0.786 4140.00  4000.000 6730.0
2      2 0.428000 -1.130          0.852 1570.00  1400.000 3770.0
3      3 0.673000 -1.540          0.882  729.00   579.000 2410.0
4      4 0.737000 -1.720          0.891  383.00   268.000 1670.0
5      5 0.745000 -1.830          0.886  220.00   134.000 1210.0
6      6 0.704000 -1.990          0.860  134.00    71.700  909.0
7      7 0.737000 -1.980          0.890   85.90    40.300  701.0
8      8 0.742000 -2.020          0.903   57.30    23.500  551.0
9      9 0.733000 -2.090          0.917   39.40    14.200  441.0
10    10 0.750000 -2.080          0.934   27.80     8.810  357.0
11    12 0.770000 -2.080          0.952   14.80     3.670  243.0
12    14 0.775000 -2.090          0.951    8.43     1.660  171.0
13    16 0.397000 -2.820          0.459    5.06     0.801  124.0
14    18 0.409000 -2.770          0.474    3.18     0.406   91.4
15    20 0.421000 -2.720          0.490    2.06     0.215   68.8

# Plot the results
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;

# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n", main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");

Is this the right way? Which soft threshold power should I select?

Any help is appreciated. thanq.

Don't worry about it too much. I would go with power 8 based on my general experience (also reflected in WGCNA FAQ) and on the mean connectivity around 50 and median around 20, which seems reasonable to me.


How to transform a DNAStringSet from the Bioconductor package Biostrings to a data frame? I am working on Mac OS X. I am using R version 3.5.1.

I have important a FASTA file into R using Biostrings::readDNAStringSet. This creates a DNAStringSet object in R.

I would like to transform it into a data frame so that I can then plot the widths with ggplot2.

as.data.frame on my DNAStringSet object returns a data frame with only one column. Is there a function that returns a data frame with columns width, seq and names present in the DNAStringSet object?

According to the documentation (?Biostrings::DNAStringSet):

  width(x): A vector of non-negative integers containing the number of
            letters for each element in x.  Note that width(x) is
            also defined for a character vector with no NAs and is
            equivalent to nchar(x, type="bytes").

  names(x): NULL or a character vector of the same length as x
            containing a short user-provided description or comment for
            each element in x.  These are the only data in an
            XStringSet object that can safely be changed by the user. All
            the other data are immutable!  As a general recommendation,
            the user should never try to modify an object by accessing
            its slots directly.

  as.character(x, use.names=TRUE): Converts x to a character vector
            of the same length as x.  The use.names argument controls
            whether or not names(x) should be propagated to the names
            of the returned vector.

This should be enough for you to create a data.frame:

dss2df <- function(dss) data.frame(width=width(dss), seq=as.character(dss), names=names(dss))

Also, if you have a huge DNAStringSet, perhaps you do not want to convert all the sequences to characters; storing sequences as characters supposedly is less efficient. For plotting it might be sufficient to have just widths.

qplot(width(dss), geom='histogram'))


Ultimate reproducibility in R? I'm looking for a convenient and reliable way to make an R analysis reproducible, either at different times or across collaborators.

Listing the package versions or a sessionInfo() output is not very helpful, as it leaves the work of re-creating the environment entirely to you.

As I see, in order to make analysis reproducible, the following has to be the same:

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

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

You're right about trying to get the R environments to be as similar as possible, as easily as possible. It's not always necessary to have the exact environment to get the same results, but it definitely helps if you find yourself struggling to replicate results.

General tips for ensuring reproducibility, not specific to R

See Sandve et al., PLoS Comp Bio, 2013. They mention a number of things that you do, like recording seed values and session info. I'll highlight a few points:

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

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

Data: Share raw and/or processed data where possible

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

Code: Comment and document your code

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

Code: Use version control

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

No manual manipulation of data

If something is done to your data, write it down

R-specific tips

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

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

Make unit tests for your work with the testthat package

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

Consider Docker images for your environment

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


Python API for working with ENSEMBL genomes Dear bioinformaticians,

I could find two python APIs to mine ENSEMBL genome databases. There is a pyensembl which is more popular and then there is genomepy which is less popular. Following is my experience with both.

1. pyensembl:

pyensembl is a fast and powerful API to mine ENSEMBL genome databases.
I have been using it for a while and in my experience it works very well. The name of the package may make you believe that it is developed/maintained by ENSEMBL. That is not true. Still, being an open source package with a realistic potential to be the default API for mining ENSEMBL genome databases.

However, there are some issues related to it's development that users may find disconcerting.

For example, documentation. I could only find a comprehensive api documentation for release 0.8.10. Current release is 1.1. Not even release notes are available. Therefore users have to rely on few lines of documentation that is available in the README.md.

Also, as far as I know, it is not peer reviewed even though it is an excellent analysis resource.

2. genomepy

genomepy loads all of the genome files into the RAM, making it immediately unusable, especially if you are working with larger genome such as human genome.

So my question is: Is there any other better python API for working with ENSEMBL genomes?

Now there's another option, the ensembl_rest module, a thin wrapper around the Ensembl REST API to simplify its usage and make it more pythonic. You can find the documentation here.

To clarify things, I'm the creator and maintainer, but still think it's a legitimate alternative.


Counting repeated kmers sequences that match at least x % of reads sequence Working on a fastQ file, I would like to get the occurrences of repeated sequences for all possible kmers of a given length that cover at least 90% of the read's length for the whole data set.

example :

  for a length 6 with the kmer "ATTGGG" and a data set  containing reads of 300 bases each

I want to obtain the occurrences of the following kmers :

ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGG......  (90 % of the read's length without  gaps)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGA.....  (91%)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGAT....  (92%)

this for all possible combinations of kmers of length = 6.

I have done some research and tested many softwares such as (KMC2,Jellyfish,Khmer,Dsk,Scturtle) but none matched what I would like to do, because I want to obtain the occurrences of the repeated kmer without any possible gaps in-between with a certain % of read coverage.

I've done a naive algorithm that works but the optimization isn't the best obviously. Is there a tool that does the trick ?

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution.

EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k*2 or longer is counted. For example in sequence ATTGGGATTGGGATTcGGGATTGGG, the script returns 15, because the second is not long enough. Counting length smaller than k*2 will be slower and more complex, but is often not the right thing to do. For example, in ATTGGGATTGGGcccccccgaaatcgatagcatcgaGGGATTGcgatc, we would not count the second GGGATTG as a tandem repeat.

On performance, for an input string of length $l$, the time complexity of this script is $O(l)$ associated with a small constant. It will be efficient and adequate to practical uses.

If you use TRF, you should run it on your sequences and then test if any tandem repeats it finds match GGGAAT – the procedure will be similar to my script, except that TRF will also find impure repeats.

// given the repeat unit length _k_, find all tandem repeats of 2*k or longer
function trf_k(k, str)
{
    var streak = 0, a = [];
    for (var i = k; i <= str.length; ++i) {
        if (i < str.length && str[i] == str[i - k]) {
            ++streak;
        } else {
            if (streak >= k)
                a.push([i - streak - k, i]);
            streak = 0;
        }
    }
    return a;
}
// find the sum of tandem repeat length with _kmer_ being the repeat unit
function trf_kmer(kmer, str)
{
    var a = trf_k(kmer.length, str), sum = 0;
    for (var i = 0; i < a.length; ++i)
        if (str.substring(a[i][0], a[i][1]).indexOf(kmer) >= 0)
            sum += a[i][1] - a[i][0];
    return sum;
}

var str = 'xATTGGGATTGGGATTGxGATTGGGATTGGGATTGGGATTGxx';
console.log(trf_kmer('GGGATT', str));


How to convert BAM file to bigWig in Python? I am manipulating some BAM files using pysam package which seems to be very fast and handy. However, I ran into problem when I am trying to generate bigWig files from BAM files in Python. I used bamCoverage from deeptools to do this previously, but it is a standalone tool that cannot be used from within my Python scripts.

Therefore, I am wondering what Python package would allow me to read in BAM files and convert them to bigWig files.

Thanks

deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:

from deeptools import writeBedGraph
from deeptools.getScaleFactor import get_scale_factor

wr = writeBedGraph.WriteBedGraph(options and input files)
func_args = {'scaleFactor': get_scale_factor(...options...)}
wr.run(writeBedGraph.scaleCoverage, func_args, output_file_name, ...)


How to align output of grep --color=always? (To QC fasta/fastq files) Grepping out short sequences from a fasta or fastq file is a really useful way to look at sequencing data. Using the option --color=always makes this even more useful, as you can visualize where the sequences appear in sequencing reads. For example:

zgrep 'GATCGATC' --color=always My_Hiseq_Data.fq.gz

To make this more useful, how can I align the output of grep --color=always such that the beginning of the highlighted parts are aligned?

Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match.

zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC

Here are some explanations:

zcat will just decompress a file.
grep will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.
first sed command will add 3 symbols before all the matches.
second sed command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of grep -Eo '.{3}GATCGATC.{0,3}'.
Finally, now that all matches are aligned, grep will color your sequence.

This code is possible to adapt to include more symbols of context.

I would recommend doing this in another language, though. E.g., in Python:

#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'

filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')

UPD

A more advanced Python version inspired by this answer with support for multiple files, gzip, etc.

#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored

def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')

# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')

if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)