Smith-Waterman traceback with large gaps I am currently trying to understand how the traceback algorithm is supposed to work for the smith-waterman algorithm as my current understanding breaks down in case of a large alignment gap.
Assume the sequences
3GA12CA
GA56CA

The obvious best alignment would be (or some other variation thereof)
3GA12-CA
-GA--6CA

However, assuming an affine gap penalty function of x => x + 2 and a binary similarity score (i.e., if equal then 1 else -1) would lead to a scoring matrix
    3 G A 1 2 C A
  0 0 0 0 0 0 0 0
G 0 0 1 0 0 0 0 1
A 0 0 0 2 0 0 0 0
6 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 1 0
A 0 0 0 1 0 0 0 2

When doing the traceback I would start with the largest value towards the lower right corner of the matrix

(6,6)
(5,5)
Now I do not know how to find the next best alignment (in this example at (3,2)) as neither going up nor going to the left nor going diagonal is going to lead to a score greater 0 meaning I have no next best value to go on.

How is the traceback function supposed to continue here?
Additionally, I read that the traceback should start at the highest score in the matrix itself. However, if we would extend the first subsequence match length by one, this would mean we miss the right alignment option completely.

This question is based on a misunderstanding. The Smith-Waterman algorithm calculates the best local segment. You have to start at the $max(H_{i,j})$ point and run the traceback as described. It will (most likely) return only a local segment that, however, has the maximum value in terms of alignment score.
To get a global alignment, another algorithm has to be used. Needleman-Wunsch would be an example.
The difference between local and global is that a local alignment only alignes a subsequence of the provided sequences whereas a global alignment gives the best alignment of the whole sequences.


how to solve “dim(X) must have a positive length” at running ComBat function in R I used ComBat() for batch effect correction in my expression data. basically, that function inputs are expression data, Batch covariate, and Model matrix for the outcome of interest and other covariates besides batch. So, I prepared all inputs based on the instruction of ComBat() function in 'sva' package. "Pheno_LMS" is my clinical data inclus batchId of my samples and "TCGA_expr_log" is my expression data. my code is below:

Pheno_LMS structure:

Sample        batchId   age  . . .
GSM71019.CEL    396     63
GSM71020.CEL    396     58
GSM71021.CEL    410     85
GSM71020.CEL    411     58
GSM71021.CEL    410     40
.
.
.
dim(Pheno_LMS)
[1] 74 37

TCGA_expr_log structure:

                   GSM71019.CEL 1      GSM71019.CEL 2        GSM71019.CEL 3  . . .
    Gene id
  ENSG00000000003     10.115170           8.628044                  8.779235
  ENSG00000000005     5.345168            5.063598                   5.113116
  ENSG00000000009     6.348024            6.663625                   6.465892
    .
    .
    .
dim(TCGA_expr_log)
[1] 42817    74

batch = Pheno_LMS$batchId
modcombat = model.matrix(~1, data=Pheno_LMS)
TCGA_expr_Co <-  ComBat(as.matrix(TCGA_expr_log),batch = batch,mod = modcombat, par.prior = TRUE,mean.only = TRUE)

but when I run that function, I got below error:

Error in apply(dat[, batch == batch_level], 1, function(x) { :
  dim(X) must have a positive length

Now, I need anybody comments about my problem. thanks

P.S: when i call batch I see:

> batch
 [1] 13 13 13 14 14 14 14 14 14 14 14 6  6  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
[29] 7  8  7  8  11 11 11 14 12 1  6  4  2  3  3  3  3  6  6  1  3  9  8  6  11 7  3  5
[57] 4  5  5  4  4  11 8  8  12 12 12 12 12 10 10 10 11 13
Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14

batch = Pheno_LMS$batchId should be batch = Pheno_LMS$BatchId, if you fix that then at least the code you showed works just fine.


How to annotate optimally a fungal genome without RNA-seq evidence? Genome information：

~50M nt
2300+ contigs
No pre-trained parameters in Augustus
There are several well annotated RefSeq genomes of this genus.
There are several RNA-seq data of this genus, but not in this species.

My current strategy:
EST evidence:

I selected many RNA-seq data that include as many species as
possible, and de novo assembled these data using Trinity.
Only the longest transcripts were saved.
All transcript files were catenated into a Total_transcript.fasta file, and cd-hit-est was used to reduce data size with options -c 0.8 -n 5 -s 0.8.

Protein evidence

All fungal proteins from SwissProt was download and reduced using cd-hit with options -c 0.8 -n 5 -s 0.8.

Proteinsets from RefSeq genomes of this genus were download, catenated, reduced using same methods.

I planed to annotate genome using Maker2 with above Evidence.
I have a doubt about my strategy: Maker2 calls Augustus to ab-initio annotation, but no pre-trained model. How should I do to solve this?
Anyone else could share your valuable suggestions? Thanks anyway.

The Maker documentation does include some instructions for training ab initio gene predictors, but it assumes an abundant EST database is available. (Assumptions about what kinds of sequences will be available for a draft genome assembly have changed drastically in the last decade.)
It may be worth exploring an iterative approach in any case. You can do a first-pass annotation with Maker using the evidence you discussed, along with ab initio predictions from Augustus (and I'd recommend SNAP as well). Each gene model will be scored by Maker with the annotation edit distance (AED), representing how well the gene model agrees with the evidence (0.0 is best). From this first pass, if you can identify a few hundred reliable gene models with good AED scores (or maybe even as many as one or two thousand), that should be plenty of training data. You could use these gene models to train a new model for Augustus, using Maker's instructions or Augustus' own instruction. (In my experience, it was quite a bit of work in either case.)


Upgrading guppy basecaller on Oxford Nanopore Gridion via PPA I am trying to update the guppy_basecaller (current version installed is 3.2.10+aabd4ec)  on a ONT Gridion device. I am initially trying to update it via the nanopore PPA, as described here (requires login access to the nanopore community)
After I issued sudo apt update, I got a GPG error, saying that the public key was not found in the system.
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://europe.oxfordnanoportal.com/apt xenial-stable InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: An error occurred during the signature verification. The repository is not updated and the previous index files will be used. GPG error: https://mirror.oxfordnanoportal.com/apt xenial-stable-gridion InRelease: The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch https://mirror.oxfordnanoportal.com/apt/dists/xenial-stable-gridion/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Failed to fetch http://mirror.oxfordnanoportal.com/apt/dists/xenial-stable/InRelease  The following signatures couldn't be verified because the public key is not available: NO_PUBKEY AD9D09AF2EBBA5A8
W: Some index files failed to download. They have been ignored, or old ones used instead

I came across a possible solution shown here, so I tried to add the key to the system with the following command:
sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys AD9D09AF2EBBA5A8

The terminal was busy for too long and, after some time I got an timeout error:
Executing: /tmp/tmp.EWnOg4cPJk/gpg.1.sh --keyserver
keyserver.ubuntu.com
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpg: keyserver timed out
gpg: keyserver receive failed: keyserver error

I then came across this post where one potential solution was to change the port, so I issued:
sudo apt-key adv --keyserver hkp://keyserver.ubuntu.com:80 --recv-keys AD9D09AF2EBBA5A8

But then I got:
Executing: /tmp/tmp.sfqFuJpVb2/gpg.1.sh --keyserver
hkp://keyserver.ubuntu.com:80
--recv-keys
AD9D09AF2EBBA5A8
gpg: requesting key 2EBBA5A8 from hkp server keyserver.ubuntu.com
gpgkeys: key AD9D09AF2EBBA5A8 not found on keyserver
gpg: no valid OpenPGP data found.
gpg: Total number processed: 0
gpg: keyserver communications error: key not found
gpg: keyserver communications error: bad public key
gpg: keyserver receive failed: bad public key

Has anyone managed to update the guppy_basecaller via PPA, as described in the nanopore community page?
An alternative would be to install via the .tar.gz file, though I am first trying to upgrade the current version from the PPA to avoid conflicts with the current installation.

This is an ONT problem, not an ubuntu problem. Keys can expire over time, and need to be updated.
The keys for the Guppy basecaller are not present in the ubuntu repository; you need to get them from ONT. This is indicated by the installation instructions for Guppy (which you've already linked to):
sudo apt-get update
sudo apt-get install wget lsb-release
export PLATFORM=$(lsb_release -cs)
###
wget -O- https://mirror.oxfordnanoportal.com/apt/ont-repo.pub | sudo apt-key add -
###
echo "deb http://mirror.oxfordnanoportal.com/apt ${PLATFORM}-stable non-free" | sudo tee /etc/apt/sources.list.d/nanoporetech.sources.list
sudo apt-get update


How to sort multiple FASTA files based on their content? So I have around 10,000 FASTA files of Influenza A virus. These files contains sequences of each of the 8 segments of the viral genome and I want to separate these files into different locations based on the content of these FASTA files.
In each FASTA file, first line has the segment number.
Eg.
>KM368312.1 Influenza A virus (A/swine/Shandong/01/2009(H1N1)) segment 3 polymerase PA (PA) and PA-X protein (PA-X) genes, complete cds

Okay imagine you have 10,000 text files in one folder and there are 8 types of these text files based on the content of them. Eg. file1.txt has segment 1 written in it, file2.txt has segment 2 written in it upto segment 8. So I want to group the files containing segment 1 together in one folder, segment 2 in other folder together and so on.

That's quite easy with cmd:
@echo off
for /l %%i in (1,1,8) do (
  for /f "delims=" %%f in ('findstr  /mc:" segment %%i " *.txt') do (
    ECHO move "%%f" "segment %%i\"
  )
)

findstr /mc:"string" *.txt searches each matching file for the string and outputs the filename, if the string was found. The for /f loop around catches those filenames, so you are able to move them. The for /l loop just iterates from 1 to 8.
See for /? and findstr /? for more information about those commands.
You probably want to adapt the filemask *.txt to your needs.
Note: I "disarmed" the move command for troubleshooting. Remove the ECHO when you are sure it does what you want.


Keep Format and Individual fields when annotating VCF with VEP I'm currently updating my Variant Calling Pipeline by switching the VCF annotating software from Annovar to VEP for a variety of regions, not least how easy it is to annotate with HGVS notation and keep datasets up to date in VEP.

For the most part everything is running smoothly, with the exception that some of the data in the VCFs is lost during annotation (and conversion to tsv). The VCFs are created with GATK's UnifiedGenotyper and include a 'Format' column where each value is 'GT:AD:DP:GQ:PL' and a column named after the Individual, which contains semicolon-separated data that corresponds to the Format column (i.e. Genotype;Allele Depth;Depth;Genotype Quality;Phred-likelihood). When I annotate with VEP none of this data is carried over to the output file as it would be in Annovar, leaving me with an annotated file that has no information on read depth, genotype or any of the other data in the two lost columns.

I've included the command I'm currently using for annotation:

./vep -i RM0108.vcf --cache --force_overwrite --tab --merged --variant_class --sift b --polyphen b --hgvs --symbol --canonical --check_existing --af_1kg --af_gnomad --humdiv --pick -o RM0108.tsv

I can't find information on this in the VEP documentation or elsewhere online. I could write something to take the relevant information from the VCF and add it to the tsv after VEP has finished running, but it seems like there may be an easier solution that I'm missing, so any help would be appreciated.

I've also posted this question on BioStars Link Here

The VEP TSV format can only keep in its own specified columns. If you want to maintain the data from your original input, get your output in VCF. It will add the VEP annotation to the INFO column, and keep all the stuff you already have there.


I want to write python program to cut a DNA sequence at an EcoRI restriction site and print the two fragments after cutting Seq =
AAAAAAAAAAATTTTTTTTTTTGAATTCCCCCCCCCCCGGGGGGGGGGG
Restriction site = GAATTC
I want to cut the sequence at GA/ATTC
I tried split method in python but It does not cut at GA/ATTC.
Please help

import Bio
from Bio.Seq import Seq
my_seq = Seq("aaaaaaaaaaatttttttttttgaattccccccccccggggggggggg")
from Bio.Restriction import *
print(EcoRI.search(my_seq))
print(EcoRI.catalyse(my_seq))

Actully the EcoRI cuts the DNA after G in restriction site GAATTC.


Get a single number representing the contact probability between a pair of genomic loci using HiC data Hello to the experts in analyzing chromosome structure data in HiC format,

I have a very basic question. I have a HiC file (specifics are mentioned below), using which I would like to obtain contact probability between a pair of genomic loci.
I obtained contact probabilities between the two loci at different resolutions available in the HiC file. However, rather than having a vector of contact probabilities, I would like to have a single number that appropriately represents the contact probability between the loci.

My guess is that there could be two ways to get a single number as the contact probability between two given genomic loci.

Should I take the highest contact probability?
or,
Depending on the lengths of the genomics loci of interest, may be the resolution that is close to that could be selected as a single number that represents the contact probability between the loci. Eg. with this strategy, if genomic loci of interest are ~5kb in length the contact probability at 5kb resolution is selected as a single number. Then, the issue is that the length of the genomic loci of interest could be lower than the lowest length of resolution available in the HiC file.
or,
Contact probabilities at different resolutions could be aggregated. Then I wonder what would be the appropriate way to aggregate contact probabilities, e.g. by multiplication or median or something else.

Just to reiterate, above these are my guess and I am not sure whats the best way.
So basically my question is: how can I get a single number that would represent the contact probability between a pair of genomic loci using HiC data?

Some specifics that might be useful in answering my question.

source of the HiC file: https://ftp.ncbi.nlm.nih.gov/geo/series/GSE63nnn/GSE63525/suppl/GSE63525_K562_combined_30.hic
organism: human (K562 cell line)
Resolutions in HiC data (bp): [2500000, 1000000, 500000, 250000, 100000, 50000, 25000, 10000, 5000]
Type of genomic loci: upstream non-coding regions of the genes
lengths of genomic loci of interest: average ~1000bp (range: 500-5000bp).
software I am using to "extract" data from HiC file: https://github.com/aidenlab/straw

The fundamental unit of a Hi-C contact is the nearest restriction enzyme cut site.
This is where each mate in a read pair from a Hi-C library comes from and is mapped to when converting sequencing data into a contact matrix.

Because the distance to the nearest cut site is variable at different genomic positions, and depth of coverage of any single locus is typically low, contacts are aggregated into bins and "normalized" to account for differences in cut site densities within those bins.
This normalization transforms the raw contact matrix into something resembling a stochastic transition matrix, where you can think of the elements of that matrix as "what is the probability that at any moment in time bin $i$ is in contact with bin $j$".

So the probability you're looking for is the $(i, j)$th element of the normalized contact matrix, where bins $i$ and $j$ contain your regions of interest.
If your regions straddle the bin boundaries, you may want to average those bins.

You get to pick the resolution of the contact matrix, but if you select a resolutiong that's too small, you run the risk of making your matrix too sparse and having no contact between those loci.
If you pick one that's too large, your matrix elements aren't really informative for your question at hand.

Given your 500-5000 bp range and median of 1 kbp, it may make sense to pick a resolution between 1 kbp and 5 kbp and use that.


Differential expression analysis for a subset of transcripts I want to check the differential expression of a specific class of transcripts (say, long non-coding RNAs) using DESeq2. Now, I know that the normalization step takes into account the total number of reads per sample. If I only input a subset of all the transcripts, the sequencing depth will change.
So my question is: can I simply input the subset of interest? Or should I input all transcripts, and subset the results later?

It's safer to import everything.  You'll want all the data for normalization and dispersion estimates.


How do I computationally predict the binding position of a DNA binding protein on the promoters of differentially expressed genes? I have a list of human genes that have been proven to be up-regulated in a disease condition of interest, using microarray analysis. I also have a protein which according to literature is said to be a DNA binding protein. I wish to be able to computationally predict if the protein mentioned binds to any of the promoter sequences of my differentially expressed genes and if yes , then at which position?

Are there any Bioinformatic tools available for this task?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you binding sites: genomic intervals where those TFs bind.

https://bioinformatics.stackexchange.com/a/2491/776

You can then map or intersect promoter regions with those TF binding sites with BEDOPS bedmap or bedops, respectively:

https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html

The use of bedmap is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.


Which tools are available to call somatic mutations on non-coding regions from whole genome sequencing data? I have a set of whole genome sequencing samples, some of which have matched normal while some not. I want to call somatic mutations on non-coding regions.

I was looking for GATK best practices (like this) but it seems it is not yet ready for production. Which alternative tools do you suggest?

Those will have to be treated as separate analysis.  You can use Mutect or VarScan for the matched samples and the tools listed in the other question for the unmatched samples however that means that you won't be able to compare the results easily.  If the inter-sample comparison is important you will have to run all samples as unmatched.  Personally I have used LoFreq for somatic variant detection with some success.  It is very sensitive so the false positive rate can be high depending on how noisy your data is.


Score predefined ChIP-seq peaks with MACS2 or equivalent I have performed peak calling on a number of separate ChIP-seq experiments and would like to harmonise the peaks from each of the experiments in order to convert my data into a matrix for further analysis. I would like to do this by taking the union of the called peaks (available in BED / narrowPeak / MACS2 format) and retrieving the enrichment score / q-value for all of the experiments on each peak. Is this possible without writing my own extension to MACS2?

You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like:

macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7

You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise.

But you can then use these "non-significant" values in combination with a bedtools map command, to map them to your predefined union of "actual" peaks.


Tool for rna/lna melting temperature prediction Is there any tool available to predict the melting temperature for rna/lna oligos with its perfect complement ?

I know this amazing website from qiagen/geneglobe that estimates the tm value.

But is there any command line tool that does this?

I can estimate the normal tm value with RNAcofold and the -T option but there is no lna option.

I searched for a while but i can only find online services.

There are several programs for simulating/predicting DNA melting but you could try them for RNA.
You can look here http://bioinformatics.org/meltsim/wiki/Main/HomePage.


Running Snakemake in one single conda env I have been experimenting a lot lately with Snakemake, I love it. Recently I also switched to using conda (--use-conda) in the way that is advertised. However, I have some issues with it, mostly related to the way we work. I work in a lab that processes data from our sequencer but also external datasets, moreover, multiple people use my pipeline.

So, with respect to conda this means: currently I have 5 envs, averaging 1.5 GB per env, this takes about 20 min to set up and it's quite some data to pull down. Especially since it is pulled down every time someone analyses data (sometime multiple times a day). This feels like quite a waste of resources. Also, I fix some package versions in the yaml files, but certainly not all (this may lead to issues when repos get reshuffled I read), so I also can't guarantee the exact environment long term (say over years, dependencies may change/update).

So I'm considering switching to a central conda env which I version and of which I export lists of all packages installed for reproducibility reasons. This will then contain all tools.

My question is, is there a nice way to do this? I had issues with putting a "source activate" line in every shell command (and how to do this for Python code under "run" anyway) before, see this: https://stackoverflow.com/questions/54024862/snakemake-conda-env-parameter-is-not-taken-from-config-yaml-file.

So, is there a way, also when using a qsub cluster, to always use a specific conda env, for everything you do?

A good practice is to provide yml file along with your Snakefile. This file should contain all the enviroment definitions that are needed to run your pipeline. The minimum content of this file would be to install snakemake itself. Any other program used in your pipeline should be included as well. For reproducibility version numbers must be used.

So, if someone would like to use your pipeline, (s)he first creates a new enviroment based on your definition:

$ conda env create -f my_pipeline.yml

Afterwards activate this environment and run the pipeline:

$ conda activate my_pipeline
$ snakemake

The conda directive in your Snakefile is only neccessary if you need to use program which dependencies cannot be resolved within the main environment. If you start snakemake with the --use-conda parameter, snakemake will create a new environment based on the provided definition file and activate it for the rule.

  Especially since it is pulled down every time someone analyses data (sometime multiple times a day).

The environment creation is only done once. If you run the pipeline a second time snakemake can activate the environment diretcly without downloading and installing something again.

The tutorial I wrote some time ago on biostars about snakemake and conda might be useful for you.


How do I Investigate and test the hypothesis on the effect of street-light regime on insects? I need advice on a suitable analysis of a dataset that has an akward design because we investigate a real-world situation. We  have two sites (A, B) where we can conduct the research. The hypothesis is: under-dimmed light, we catch less number of insects than under-fully lit street lights.

The two sites each have ten street lights that were sampled during 16 nights of dimmed light and during 16 nights of full light. During these total of 32 nights, we sampled insects at each street-light pole (10 per site, 20 in total). For A and also at B, we also have temperature as a variable (nightly average, max, min) - a measure per site not per light pole.

I would like to fit a model: number of  insects ~ dimming, temp_mean, temp_min, temp_max, site

I  include "site" as categorical variable and would like to use the 32 nights (16 dimmed, 16 fully lit) as 32 unique samples (each night as one sample). I am not sure whether this is "statistically sound". If not, does anybody have an idea on how to evaluate such data in order to prove the hypothesis? Would you go for t-test, PCA or something else?

streetlight <- data.frame(site=c(rep("A",16),rep("B",16),rep("A",16),rep("B",16)),
                          dimming=c(rep(1,32),rep(0,32)),temp_mean=runif(32, min=20, max=24),t
                          emp_min=runif(32, min=12, max=19),temp_max=runif(32, min=25, max=28),
                          insects=round(runif(32, min=5, max=20),0))

The hypothesis can be examined by t-test. You can group the data (number of catches) into  dimmed light, and fully lit street lights and proceed with t-test to assess whether dimmed light results  in less en-catchment or higher en-catchment or equivalent encatchment.  Alternatively, a linear model can be applied subject to certain conditions  here. PCA is not  appropriate here.


What is the perferred method of optimization or energy minimization of small molecules downloaded from PubChem? I will be docking a small library of molecules that I have downloaded from PubChem using AutoDock. I thought why not minimize their energy before docking using ChemDraw. I am thinking; what is the best method of optimization/ energy minimization of small molecules. MM2, MMF94 or other? And why?

The difference between different force-fields is not going to be major, it is the side steps which are.

When

If you are starting from a SMILES string, optimisation is a must obviously.

If you are using a 3D conformer from PubChem or even an actual sub-1 Å crystal structure from CSD, optimisation is nice for consistency.

Which

MMF94 is a solid choice. RDkit offers MMF94 (Merck force field 1994) and UFF (universal force field). The latter has some extra features, such as not kicking the bucket at atoms that aren't your usual suspects. ChemDraw offers MM2 and MMF94 IIRC, the former is very old. g16 in Gaussian (paid software) uses Gaussian Electrostatic Model, which is polarisable and will solve partial charges in a nice way (which I'll get to in a minute). OpenBabel uses MMF94 under the hood. You can even optimise with macromolecular forcefields, such as AMBER and CHARMM.

Given that the AutoDock uses a flexible ligand it should not matter too much anyway.

What else?

Autodock does not require conformers by virtue of the flexible ligand business (the sampler does rototorsions). Other programs might if they use rigid ligands (e.g. Fred) or they use a mixed approach and use a conformer library for speed (e.g. Rosetta ligand_dock). Conformer diversity is actually a bigger issue as conformers that are only slightly contorted, often even less than kB (0.6 kcal/mol), may not be sampled.

The thing that messes things is incorrect protonation. OpenBabel can protonate at a given pH, RDKit cannot. Not an issue for ligands with drug screens (cf. Lipinski rule), but is for natural substrates. However, for the protein it is essential to make sure acidic and basic residue are correctly protonated and that histidine have the correct tautomer —there are tools, c.f. MOE, that can do it but by far the best approach is reading up on the mechanism.

The next thing is dodgy Gasteiger Marsili partial charges —for Autodock 4 (and most other full/hybrid physics-based docking programs), partial charges are important: after all PDBQT has a Q in it for charge! For prepare_ligand4.py you give it a mol2 file. A major consideration is what format you are using for your molecules. This sounds a very silly thing to point out, but you'd be surprised how many trip up. An mol file (or an sdf file, which is the same but multi-entry) does not have an official partial charge column (although they are often added as comments after M END). So whatever force-field you chose make sure that the program you use gave you partial charges are not formal charges. This does not quite apply to AutoDock Vina, which is charge independent and more reliant on machine learning derived factors.

Verdict

Open babel is nice in that it can start with a SMILES, handles multiple formats including mol2 and PDBQT, can protonate at pH 7, uses MMF94 and does Gasteiger Marsili charges (and its confab algorithm is very good, not need here though).


Canu assembly not making a single consensus? I've downloaded reads from this BioProject. Using canu with default parameters (no correction), I've got 4 contigs, none of which really look like the reference plasmid here.

The command I used was:

canu -p ip40a -d ip40a_assembly -useGrid=false -maxThreads=6 -genomeSize=175k -nanopore-raw reads.fastq

Granted, the original publication for this sequence used minimap/miniasm but my questions are:

A. Why do I get multiple contigs? Doesn't canu generate a single assembled sequence?

B. This is a plasmid but canu says suggestedCircular=false. What's going on here?

Looking that the link you provided only 30% of the reads appear to be e. coli.  I would suggest filtering the reads for those that align to e. coli and see where that gets you.  Did the original publication describe the method they used for data cleaning?

A:

https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR3385327

Unidentified reads: 68.82%
Identified reads: 31.18%
cellular organisms: 30%
Bacteria: 30%
Proteobacteria: 28.09%
Gammaproteobacteria: 19.7%
Enterobacterales: 13.3%
Enterobacteriaceae: 12.54%
Vibrionales: 0.01%
Viruses: 1.17%

B:

Without a good assembly canu won't be able to determine if the assembled genome is circular.  I would ignore this until you have a good example.


Find most similar genes to a set of genes of interest in time series RNA-seq data I have gene expression data (RNA-seq) for 30 different time points (from 0 to 60 min each 2 min).

I have a set of 8 genes that behave similarly (although not identically) and I want to find the top X genes that behave most similarly to this set of genes in terms of their time-series expression profile.

I have tried so far to do unsupervised clustering (mfuzz) and check for genes that lie in the same cluster as my genes of interest. However, this approach relies on specifying a predefined number of clusters and since these genes are not exactly equal, sometimes they don't fall in the same cluster, depending on initial parameters.

Because of this, I have thought about doing it the other way around and look for genes similar to my set of genes. However, I haven't found any method/resource for doing that so far.

So, is there a way to find genes with a time-series profile similar to a predefined set of genes?

You can use GSVA to summarize the information of your 8 genes and look for the time-serie change using limma or DESeq2.

If you have enough samples for each time point you could also do a weighted gene co-expression network analysis (WGCNA), and see if those clusters in each time point are co-expressed similarly (But it can be tricky and you will need to define a lower bound of the clusters to 8 or so). You could also divide your time series in two point at the beginning and at the end, and see if there are differences or not.

If you want to get a list of genes similar to your set of genes you can use GOSemSim or BioCor (Disclaimer I am a contributor to the first, and author of the second), which can calculate similarities between genes based on the gene ontology and on pathways respectively. For this approach you'll need to calculate the similarity between all the genes and then select those that are closer to your genes of interest.

I'm not aware of any database/website with time-series profiles of genes (I assume that if there would be any would be on humans or some simple model organism).


Combine multiple species trees I have five species trees run using starBeast and Beauti. Here is a link to the paper discussing both,
http://journals.plos.org/ploscompbiol/article/file?type=supplementary&id=info:doi/10.1371/journal.pcbi.1003537.s006
 Each of these trees has the same individuals, and each locus used (10 loci per tree) has individual models set to it. All the trees used a relaxed log clock normal. I would like to be able to combine these trees into one. I have attempted to find a way to do this within Densitree and starBeast, but to no avail. When all loci are run with all trees, the program simply dies, leading me to believe that running the trees together will all parameters is computationally impossible for my hardware. Is there a way to combine these trees?

Although I've never tried it myself, I think what you need can be achieved through the LogCombiner program from the BEAST package -- assuming the species trees files are comparable (ran with same parameters, same species names, etc.) I am also assuming that you meant 'five species tree files', each file with the posterior sample of species trees.

However, I do believe you should try to run all loci within the same analysis, or at least do some prior check to see if the distributions between data sets are comparable. LogCombiner assumes independent runs from the same data.

I would also recommend taking a look (and following) the BEAST google group (which includes starBEAST discussions).


Searching for proteins by primary and secondary structure I need to find peptides/proteins with a certain sequence of amino acids in a certain type of secondary structure. For example, "ALA PHE GLY" in alpha-helix. Is there any tool that does that?

Your request is rather specific, so there is no tool...

It is a straightforwards request though. I'd impose the requirement that the structures be solved and no predicted secondary structure. You can do it two ways if using Python.

PDB

After downloading all PDB structures, in python with PyMOL (installable with conda, tricky without) you'd do something like

import pymol2, os
from multiprocessing import Pool

def find_motif(file):
    motif = 'AFG'
    try:
        with pymol2.PyMOL() as pymol
            pymol.cmd.load(file)
            seq = pymol.cmd.get_fastastr()
            if motif not in seq:
                return None
            d = {x: []}
            # str.find will cause issues if there is a gap also structure may not start at one.
            pymol.cmd.iterate('resn ALA+PHE+GLY', 'x.append((resi, resn, ss))', space=d)
            sequence = -1
            i = 0
            for residue in d['x']:
                if residue['resn'] = 'ALA':
                    sequence = 0 # lazy enum
                    i = int(residue['resi'])
                elif sequence != -1 :
                    if int(residue['resi']) == i+1 and residue['resn'] = 'PHE':
                        sequence = 1
                    elif int(residue['resi']) == i+2 and residue['resn'] = 'GLY':
                        return (file, i)
                    else:
                        sequence = -1
    finally:
        return (file, None)

if __name__ == '__main__':
    p = Pool(12)
    structure_files = [f for f in os.listdir('pdb') if os.path.splitext(f)[1] == '.pdb']
    verdicts = p.map(finder, structure_files)
    print([(f,v) for f,v in verdicts if v is not None])

The end result will be a list of PDB with where in the sequence (in that PDB code).
(I wrote this in the reply box, so there may be a typo or two).
PyMOL has a selection algebra called pepseq, which matches missing atoms, but I was a whimp and played it safe by being more verbose. Also, a for loop can be used instead of an iterate. But I think it is slower.

atoms = pymol.cmd.get_model(f"pepseq {target} and name CA")
for atom in atoms.atom:
      if atom.ss == ' ...

Uniprot

In Uniprot, if a protein is crystallised, there will be the feature entries corresponding to the secondary structure. It is sheet or helix. No turns, pi-helices or 3_10 helices. There is no Biopython method for parsing Uniprot XML. I wrote myself a parser and it is not too horrid thing to do —although I monkeypatched elementtree. Alternatively, xmlschema.XMLSchema.to_dict can be used to convert the XML to a dictionary. The advantage of Uniprot is that the residues ids will not require correcting. The disadvantage is parsing XML.

My parser is specific for my application and it is poorly documented (so would not recommend it) but but... I can share it and human data. I'd be used like:

 import os
 from michelanglo_protein import ProteinCore
 ProteinCore.settings.startup(data_folder='protein-data')
 p = ProteinCore().gload(file=os.path.join('gpickle','Q86V25.gpz'))
 #print(p.asdict())
 if 'helix' in p.features:
        n = p.sequence.find( ...
        for ... p.features['helix'] #('helix', 'turn', 'strand')
        if x < n < y ...


stacked bargraph with error bars I'm attempting to plot a stacked barplot with ggplot2 with this code:

ggplot(CC, aes(x = Condition, y = Percent, fill = Cell_Cycle))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))

Question: How do I go about adding error bars to the graph?

I tried this code:

ggplot(DF, aes(x = Condition, y = Percent, fill = Cell_Type))+
     geom_bar(stat = "identity")+
     geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5) + geom_errorbar(aes(ymin=Condition, ymax=Condition), width = .2))

But got this error:

  Error: Cannot add ggproto objects together.

output of dput(DF)

structure(list(Cell_Type = structure(c(1L, 2L, 3L, 4L, 5L, 6L,
7L, 1L, 2L, 3L, 4L, 5L, 6L, 7L), .Label = c("Fibroblast", "T cell",
"Macrophage", "Tumor", "Islets of Langerhans", "Endothelial",
"B cell"), class = "factor"), Condition = structure(c(1L, 1L,
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("treated",
"untreated"), class = "factor"), Freq = c(6051L, 1892L, 1133L,
657L, 116L, 868L, 832L, 5331L, 3757L, 1802L, 835L, 287L, 704L,
256L), Percent = c(52.3941466793662, 16.3823707680319, 9.8103731924842,
5.68880422547407, 1.00441596675037, 7.51580223395965, 7.20408693393367,
41.0962072155412, 28.9623805118717, 13.8914585260561, 6.43694110391613,
2.21245760098674, 5.42707369719396, 1.97348134443417)), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"), row.names = c(NA, -14L), groups = structure(list(
    Condition = structure(1:2, .Label = c("treated", "untreated"
    ), class = "factor"), .rows = list(1:7, 8:14)), row.names = c(NA,
-2L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE))

If you want to indicate the error or the uncertainty for each cell type, you need to replicate your experiments, ccount each cell type in each replicate, calculate the average percentage of each cell type in each condition and the associated standard deviation.

To answer your question about how to plot a stacked bargraph with error bar, first this option is not integrated into ggplot2 (see this discussion from the creator of ggplot2: https://github.com/tidyverse/ggplot2/issues/1079)

So, you need to calculate the position of each error bar by hand. You can do that by using dplyr. Here, I set the same error bar for everyone as you don't know what is the standard deviation in your data but at least you can observe the structure of the dataframe:

library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  mutate(Cell_Type = factor(Cell_Type, levels = c("B cell","Endothelial","Islets of Langerhans","Tumor","Macrophage","T cell","Fibroblast"))) %>%
  group_by(Condition) %>%
  mutate(SDPos = cumsum(Percent)) %>%
  ggplot(aes(x = Condition, y = Percent, fill = Cell_Type))+
  geom_bar(stat = "identity")+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_stack(vjust =  0.5))+
  geom_errorbar(aes(ymin = SDPos-SD, ymax = SDPos+SD), width = 0.3, position = "identity")

As you can see, it's pretty ugly... so, if you need to add error bar, I will rather advised you to do a dodged barchart as the following:

library(dplyr)
library(ggplot2)
DF %>% mutate(SD = 2) %>%
  ggplot(aes(x = Cell_Type, y = Percent, fill = Condition))+
  geom_bar(stat = "identity", position = position_dodge())+
  geom_text(aes(label = paste(round(Percent,2),"%")), position = position_dodge(0.9), vjust =-2, size = 3)+
  geom_errorbar(aes(ymin = Percent-SD, ymax = Percent+SD), width = 0.3, position = position_dodge(0.9))+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_y_continuous(limits = c(-2,60))

To my opinion, it makes easier the comparison of the distribution of each cell type between conditions.

Hope that it answers your question.


Does module size decrease as the number of samples increase in WGCNA? We are working with a lot of samples, more than 50,000 people! We are using WGCNA to get modules and it seems like there are less number of module genes the more people we use. We are following the tutorial guidelines. We set the min modulesize = 30 and deepSplit = 4 (max sensitivity) However, we are getting very small modules with numbers almost 30 for all. We want to know if using this many samples creates an artifact and the variation between the samples may reduce the correlation seen between the genes? If that is the case, do you suggest any way we can fix that? Also do you think z-scaling expression data improves the module detection?
Thank you so much.

The biggest difference a large number of samples makes is that you can usually decrease the soft thresholding power. If you use power 6 (or 12 for a signed network), try decreasing it to 3 or even 2 (6 or 4 for signed network). Apart from that, it could also be that the co-expression modules in your data are simply smaller than what people usually find, nothing apriori wrong with that.


Visualisation of long read RNA-Seq splicing I have a dataset of Oxford Nanopore cDNA reads. Many of my reads are full-length or close to full-length transcripts, and I and am interested in examining alternative splicing. For this, I would like to begin by visualising my reads and comparing variants qualitatively.

I have tried visualising my data in both IGV and SeqMonk, but neither has given a satisfying result.

IGV

IGV shows some links between aligned sections, shown as fine blue lines and thick black lines. I cannot work out what the difference between these lines is, and even more confusingly, we have numerous alignments where two aligned exons come from the same sequence but are not joined by any line. This means that, to confirm or deny a spliced variant, I need to manually examine the read ID of each exon.

SeqMonk

SeqMonk gives a much cleaner visualisation, but unfortunately shows no links between aligned exons, and worse, I cannot find the read IDs. This means there is functionally (as far as I can see) no way to tell which exons come from the same read.

EDIT: Sashimi Plot

As recommended in this answer, I have tried using IGV's Sashimi Plot. Unfortunately, I believe the high error rate and relatively low coverage is causing this to create somewhat confusing output. The Sashimi Plot shows mostly nearly every junction to have just a single read supporting it.

Is there a visualisation tool which allows simple examination of splicing of full-length transcripts?

DNASTAR's software is for purchase, but high quality. GenVision Pro does genomic visualization, including Sashimi plots.

Edit: not sure why this answer is being downvoted, unless it's because the software isn't free. OP has tried IGV and SeqMonk, I mentioned an alternative he might not have heard of.

Here is a video demonstrating the use of Sashimi plots in GenVision Pro:

http://www.dnastar.com/t-support-videos.aspx?video=YJvcERoSIsg


Handling sample identity in aggregated 10x libraries? cellranger aggr can combine multiple libraries (samples), and appends each barcode with an integer (e.g. AGACCATTGAGACTTA-1). The sample identity is not recorded in the combined matrix.mtx file.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate#gem_groups

How can we keep and assign sample information to each cell after loading the data into R? (e.g. Seurat::Read10X())

You have to keep track of the file order you used for aggr.  The suffix number represents which represents the aggregated sample. You could store this information in a text file and load it into R independently of the Read10X function and combine it with the cell names to get their sample names. From the page you link:

  Gem Groups

  This is an integer that is appended to each barcode in the
  gene-barcode matrix. For example, AGACCATTGAGACTTA-1 and
  AGACCATTGAGACTTA-2 are distinct cell barcodes from different
  libraries, despite having the same nucleotide sequence.

  The numbering of the GEM groups will reflect the order that the
  libraries were provided in the Aggregation CSV.


What kind of "gff" format does bioawk parse? I was wondering if I could use the gff parsing capability of bioawk to facilitate the parsing of gtf files, and I looked at the following help message:

$ bioawk -c help
bed:
    1:chrom 2:start 3:end 4:name 5:score 6:strand 7:thickstart 8:thickend 9:rgb 10:blockcount 11:blocksizes 12:blockstarts
sam:
    1:qname 2:flag 3:rname 4:pos 5:mapq 6:cigar 7:rnext 8:pnext 9:tlen 10:seq 11:qual
vcf:
    1:chrom 2:pos 3:id 4:ref 5:alt 6:qual 7:filter 8:info
gff:
    1:seqname 2:source 3:feature 4:start 5:end 6:score 7:filter 8:strand 9:group 10:attribute
fastx:
    1:name 2:seq 3:qual 4:comment

I see that there are 10 fields defined for gff parsing.

However, when I look at the pages for gtf and gff in the ensembl website (gff/gtf and gff3), all have 9 fields.

I'm curious about these "filter" and "group" fields. What are they meant for? Are they extracted from some of the columns mentioned in the above pages?

I would consider the description there a bug. The filter is actually the strand, strand is the frame, group is the attribute, and attribute does nothing. These are really meant to be the 9 columns.

Edit: There's a bug report related to this.

Edit 2: I've made a pull request to clarify this and fix the aforementioned bug report.

Edit 3: I realized that I never directly answered the title of your question (mea culpa). bioawk itself will work with gff, gff3, or gtf files. It really is just treating them as tab-separated files with named columns (this is surprisingly convenient, since it's a PITA to remember what column does what).

Edit 4: The PR has been merged. If you install from github then you'll see corrected field names.


How to correlate two zero inflated bedgraph-like signals? This question pertains to iCLIP, but it could just as easily be ChIP-seq or ATAC-seq or mutation frequencies.

I have iCLIP read counts across the transcriptome and I wish to know if the signals are correlated - that is, where one of them is high, is the other likely to be high.

Often when dealing with such data (e.g. iCLIP data) we know that the data is generally sparse - that is at most positions both signals are zero and this is correct, and also zero-inflated - that is many bases that "should" have a signal are missing that data. So just calculating the Spearman's correlation is likely to give an artificially low value.

What might be a way to asses the association? I should add that the aim is the assess association of binding patterns within genes, rather than (or as well as) between genes.

Things I have thought of:

Apply some sort of smoothing to the data (eg a rolling mean). Remove any bases with 0 in both samples. Compute the spearmans.
Calculate the average pairwise distance between every read in sample one and every read in sample two. Compare this to data where the reads have been randomised within genes.

In the first case removing all bases with 0 in both samples seems wrong. But if 99.99% of all bases have zero in both samples, then this seems like its necessary for Spearman.

In the second case, the result seems like it would be non-intuitive to interpret. And also calculating this would be massively computationally intensive.

Honestly I'd just use multiBigwigSummary and then plotCorrelation from deepTools for this, but I'm a bit biased. There, the idea would be to consider each gene as a unit (you could instead use bins, but I don't think that would as nicely do what you want), namely by giving the tools a BED or GTF file input. It would then calculate the average signal in each gene/transcript and you could do your spearman's correlation. Features with 0 in all samples could optionally be removed (plotCorrelation --skipZeros).

While you certainly could go the whole 9 yards and use per-base comparisons, that seems a bit overkill and I suspect that it won't really yield appreciably more information (especially when one considers the additional time overhead).


Get RefSeq accession numbers with versions Google searching for NM_002084 gives the following result:

NM_002084.4

This, I assume, is the latest version v4, hence the .4 suffix.

Searching for previous versions I get the following results, along with notes saying it was updated or removed.

NM_002084.3

  This sequence has been updated. See current version.

NM_002084.2

NM_002084.1

  Record removed. This record was replaced or removed.

Using biomaRt I can get the latest(?) version as follows:

library("biomaRt")

# define db
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

# get refseqs
getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = "NM_002084"),
      mart = ensembl)
#  refseq_mrna chromosome_name transcript_start transcript_end strand
#1   NM_002084               5        151020438      151028992      1

But querying for specific versions gives nothing:

getBM(attributes = c('refseq_mrna',
                     'chromosome_name',
                     'transcript_start',
                     'transcript_end',
                     'strand'),
      filters = c('refseq_mrna'),
      values = list(refseq_mrna = c("NM_002084.1", "NM_002084.2", "NM_002084.3", "NM_002084.4")),
      mart = ensembl)
# [1] refseq_mrna      chromosome_name  transcript_start transcript_end   strand
# <0 rows> (or 0-length row.names)

Question: How can I get all versions (preferably using R)?

I don't believe this is possible using biomaRt, nor using AnnotationHub.

I have two suggestions, neither of them very satisfactory. First, you can specify an Ensembl archive for biomaRt, for example:

mart72.hs <- useMart("ENSEMBL_MART_ENSEMBL", "hsapiens_gene_ensembl",
                      host = "jun2013.archive.ensembl.org")

Of course, that requires that you have some idea of the date for each accession version and that the archives span that date - so not especially useful.

The other option is to access EUtils using e.g. rentrez, which does allow search by version number:

library(rentrez)
es <- entrez_search("nuccore", "NM_002084.1 NM_002084.2")
es$ids

[1] "4504104" "6006000"

So knowing the accession, you could simply append 1, 2, 3... to it, run the search and see if UIDs come back, then get them using entrez_fetch.

EDIT: here's a quick and dirty function which takes an accession as query, appends version = 1, fetches the ID, then increments the version and repeats until no more results are returned. It is not well-tested!

getVersions <- function(accession) {
  require(rentrez)
  ids <- character()
  version <- 1
  repeat({
    es <- entrez_search("nuccore", paste0(accession, ".", version))
    if(length(es$ids) == 0) {
      break
    }
    ids[version] <- es$ids
    version <- version + 1
  })
  ids
}

Example:

getVersions("NM_002084")
[1] "4504104"    "6006000"    "89903006"   "1048339180"


How to convert a phylogeny to a dendrogram in R I want to convert my phylogeny into a dendrogram so I can use it with dendextend in R to produce a tanglegram. I have made some progress but I keep encountering errors, see below:

library(ape)
library(dendextend)

Tree <- rtree(10, rooted=F)

Tree <- read.tree(file="clipboard", text=NULL)
# test <- as.dendrogram(Tree)
## Error in ape::as.hclust.phylo(object) : the tree is not ultrametric

Tree.ultra <- chronos(Tree)
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not rooted

Tree.ultra$root.edge <- 0
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not binary

This thread seems to be the closest to answering my issues but the workflow is different and Im not sure how to address this binary issue. I have tried using as.hclust, as.hclust.phylo and hclust.

For your example to work you should replace the rtree function with the rcoal function, see example below.

library(ape)
library(dendextend)
Tree <- rtree(10, rooted=F)
is.ultrametric(Tree)
[1] FALSE
is.binary.tree(Tree)
[1] TRUE

Tree_rcoal <- rcoal(10, rooted=F)
is.ultrametric(Tree_rcoal)
[1] TRUE
is.binary.tree(Tree_rcoal)
[1] TRUE

Tree_rcoal2 <- rcoal(10, rooted=F)

tanglegram(Tree_rcoal, Tree_rcoal2)

To make use of a file in newick format, you can try to convert with the chronos function.

s <- "owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
is.ultrametric(tree.owls)
[1] FALSE
dendrogram <- chronos(tree.owls)

Setting initial dates...
Fitting in progress... get a first set of estimates
         Penalised log-lik = -24.04992
Optimising rates... dates... -24.04992
Optimising rates... dates... -24.04908

Done.
is.ultrametric(dendrogram)
[1] TRUE


Why are my kallisto and salmon results differing so much just for lncRNA transcripts? I am running some analysis on an RNA-seq dataset. I have a list of transcripts that are potential lncRNA for which I ran both Kallisto and Salmon aligners. The input data for index building and quantification includes mRNA as well as these potential lncRNA.

Commands for one replicate:

Salmon

salmon index --keepDuplicates -t transcripts.fasta -i index_directory/ -k 31
salmon quant --seqBias --gcBias -i index_directory/ --libType ISR --auxDir index_directory/ -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz -o output_name

Kallisto

kallisto index -i kallisto_index.idx transcripts.fasta
kallisto quant -i kallisto_index.idx -o output_name --rf-stranded sample_R1.fastq.gz sample_R2.fastq.gz

Then I used the tpm values, averaged replicates (3) and plotted the number of transcripts per condition that had a tpm > 2. Same result for tpm > 5.

The following graph shows the for mRNA:

And for lncRNA:

Again, the analysis were done with both lncRNA and mRNA transcripts together in the same fasta file.

The overlap is not shown here but is always smaller than the count for Kallisto, for the lncRNA.

In literature both tools seem to have similar results. I cannot find a reason for this. If there would be an error somewhere I would expect to see also differences in the conding RNAs which does not happen.

There has been a detailed comparison of kallisto and salmon for lncRNA transcripts by Zheng et al https://www.biorxiv.org/content/biorxiv/early/2018/01/09/241869.full.pdf

It shows them to be near identical. This suggests there is some problem with the way you made the comparison/plot.


Edit FASTA header using sed I need to rename the following headers from a FASTA file. Something like:

>ENST00000488147.1::chr1:29570-30013(-)
>ENST00000405392.5::chr2:24919186-24920186(-)

into different format like this:

>chr1:29570-30013(-)_ENST00000488147.1
>chr2:24919186-24920186(-)_ENST00000405392.5

How can I rename the FASTA headers, perhaps using sed or awk?

awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa

An explanation of the non-obvious bits:

BEGIN{FS="::"} split columns using :: as the designator
$1~">" If there's a > in the line
substr($1,2) Trim the > that was at the beginning of the line.


How to provide a GFF/GTF file to Chromonomer I am trying to use Chromonomer. It should be possible (according to the software page and a paper) to add a GFF/GTF file to the analysis.

Alas, neither the built-in help nor some light source code browsing were enough to find a way of passing the file to the software.

How can I use a GFF or a GTF file with Chromonomer?

I think the part of the code that should allow you to do this is in chromomer.cc.

Lines 2857-2983 form the function parse_command_line. Nowhere is there any reference to anything resembling parsing a GFF file, so I think the documentation is not in line with the code at this point.

For reference, I looked at version 1.07. It would be helpful if the code were on github.


Asking for Docking programme I'm tring to do my M.Sc. research and I Have to do Docking to save some money as I cannot try all the compound that I'm working on.

what is the best Protein-Ligand Docking programme?

And anyone has any tutorial for such programmes.

It's case specific
Docking small compounds well is problematic and you need to consider whether:

if a single conformer is okay
if rigid sidechains are okay
if implicit water is okay
if a non-polarisable model is okay
if a non-covalent interaction is okay
if you are okay with a score that is a rough ΔΔGibbs approximation of the interaction (see this SO answer) or do you want something more like a K_d

Do note that for any experiment you do to have any validity you need to have controls. Do you know what can and cannot bind already? It is imperative that you score these too.
Free
I am assuming that you want to use a free program. Virtual compound screens (VCS) are a big business in pharma. Therefore, a lot of really good software is not free, such as ICM-Dock, which outscores AutoDock discuss below in most tests.
AutoDock: simple cases
AutoDock Vina is easy to use and has a PyMOL plugin (I believe), but the peptide, in both its sidechains and its backbone, is rigid, which is not good. The ligand rotamers are generated for AutoDock at runtime, but there is a hardcoded low limit, it technically isn't a rigid ligand, but nearly.
It is the best-scoring free rigid peptide docking program.
There are lots of tutorials on using Autodock, such as the official one or on YouTube. As shown by this Bioinformatics Stack Exchange question Autodock 4 (different score fxn) has several steps, each important.
Rosetta: very complicated cases
I prefer the harder to use Rosetta Ligand_dock which uses conformers (if the ligand is parameterised correctly, say they were generate via Open Babel conformers) and allows sidechains to repack.
Water
Both tools mentioned use implicit waters, which works badly for highly hydrophobic ligands.
That is the waters are simulated as a homogeneous field that affects the ligand and protein. But some solutions are possible, in the case of AutoDock WaterDock variant exists —although this is not polarisable water.
For more advanced cases there are MD methods such as dynamic undocking, which are completely non-trivial. All these methods mentioned use classical force-field calculations (Amber, Charmm, Talaris etc. are different models), but in MD noise (think of it as environmental heat) is present in the calculations, so the ligand comes unbound (and is actually pulled out). You can also get better results by refining your bound ligand with QM-MD (e.g. gauss of gaussian).
Choices...
This is a lot of options. So a first question is are your ligands man-made drugs/fragments (i.e. obey Lipinski's rule of 5) or metabolites? Generally the former are rigid and very hydrophobic (greasy) (and require better force fields), while the latter are flexible and hydrophilic.
Minor considerations
Also to consider is your entranceway to your protein as the Michaelis constant is dictated in some enzymes by the substrate going in, not by its binding (e.g. P450s).
If the ligand have a shared backbone that is fixed there are some tricks you can do to cut your compute time and increase the reliability.
If the ligands are bound not all methods are possible. Rosetta ligand_dock can but it is a bit fiddly to parameterize.
RDKit skills
One thing also to consider is your programming skills. Can you use RDKit? In most cases downloading the sdf off PubChem and running OBabel or similar to generate conformers is not possible and you have be familiar with a python package called RDKit, which is really powerful but non-trivial to use.
Additionally, docking scores can be strongly augmented with ligand-only values (such as logP, TPSA, molar refractivity, molecular weight, QED, Bertz topology etc. etc) that are generated from RDKit.
NB. The scoring function in Autodock Vina (not 4 though, as that is Amber FF) partially takes into account.
Ligand_dock details
Ligand_dock is an older application of the Rosetta suite and the current philosophy is to use a specific Rosetta script, however, it works fine and is actually quicker to implement. Here are some links:

manual of ligand_dock
tutorial for ligand parameterisation
PubChem for your SDF files
online tools for small molecule manipulation
example of ligand creation for docking with RDKit


UCSC liftover cannot produce bed format output I was trying to use UCSC liftover tool to convert human coordinates to mouse. While most of the files worked, one file didn't:

chr16   86430087    86430726    Human_chr16:86430087-86430726___element_1___positive____neural_tube[12/12]___hindbrain_(rhombencephalon)[12/12]___limb[3/12]___cranial_nerve[8/12]
chr16   85620095    85621736    Human_chr16:85620095-85621736___element_2___negative_
chr16   80423343    80424652    Human_chr16:80423343-80424652___element_3___negative_
chr16   80372593    80373755    Human_chr16:80372593-80373755___element_4___positive____neural_tube[6/10]___hindbrain_(rhombencephalon)[10/10]___midbrain_(mesencephalon)[10/10]

Does anyone know why this bed file can only be converted to chr:start-end format but not bed4 format?

I figured out the issue by myself. It turned out that UCSC liftover doesn't like the character ":"

I removed ":" and it worked


Can foldx position-scan produce different result with different lists of mutations? I am running foldx position scan for GFPmut2 (PDB: 6GO9) and I use a config file with the following template:

command=PositionScan
pdb=6GO9.pdb
pdb-dir= <pdb location>
output-dir= <output directory>
output-file= <file name>
positions= <list of positions, comma separated: QA204R,KA107R...>

I had created three different config files with different position list.

Some of the positions were present in more than one files. Foldx produced different estimates of ΔΔG for these positions from these different config files. Example:

File1-LYSA41K   0
File1:LYSA41R   -0.353758
--
File2-LYSA41K   0
File2:LYSA41R   -0.0356187
--
File3-LYSA41K   0
File3:LYSA41R   -0.0031362

Another example

File1-LYSA52K   0
File1:LYSA52R   0.0057333
--
File2-LYSA52K   0
File2:LYSA52R   -0.00809613
--
File3-LYSA52K   0
File3:LYSA52R   0.0077189

I can certainly run foldx with a union of the three list of positions but I don't understand how this difference can even occur. Shouldn't each position scan be independent?

Minute values.
So just a reminder that the ∆∆G units is a kcal/mol. Water bounces around and collides with different energies following a Maxwell-Boltzmann distribution, where the average collision at 25°C is 0.59 kcal/mol (thermal energy constant).
A hydrogen bond reduced the ∆G by 1–3 kcal/mol. Other interactions such as π-π, π-S are on the 1–2 mark. So the difference you are seeing is totally negligible.

Heuristic.
The reason why you see the differences is that the force field calculations are heuristic and are not testing all perturbations —it is not only sampling a limited amount of conformers of side chains (kind of like PyMOL would be if it were FF based), but neighbour perturbations. So sometimes you hit a wee energy minima. Hence why in some protocols replicates are required —not mutation scanning mind you.

Bonus: Protein pI.
Lastly, a caveat about pI of protein, you probably already know this, but best be safe. Your scores are all low, so I am guessing these surface mutations? Improving thermal stability is great for increasing the T_M measured by DSC or catalysis if you are scanning a transition state as done for theozyme design, but it does not mean you get better solubility. In this case, you are introducing positive charges, which will alter the pI of the protein upwards and most protein have a pI less than 7 and the closer it is to pH 7 the more they crash out. If you are designing a high pI GFP, Rosetta supercharge does FF calculations to alter the pI of protein, which might be worth a gander.


Best data repository to publish a large 'plotly' table containing all annotations on a transcriptome? I am not sure whether this is the best SE community for this question.

I want to publish a transcriptome paper, along with interactive materials enabling readers to peruse the data behind the discussion. The R package 'plotly' enables rendering online-publishable tables that can be filtered & searched for specific entries.

I would like to publish such a table as supplementary to a paper (e.g. preprint version) however usually data repositories allow for uploading files and text.

Would anyone here know an official data repository with a format such that would enable publishing online an interactive table?

A good example of such a table is given here.

Thanks in advance.

You may try Code Ocean. Currently, it is not that easy to use IMHO, but one of their tech leads promises to substantially improve it towards the year end.


RNA_Seq data aligned used uniquely or multi mapped reads impact on result interpretation I have some transcriptomic (Whole) sequencing data that I should analyse. I would like to do raw data alignment to a reference genome taking into account the multi mapped reads and uniquely mapped reads. I would like to know the consequences on my transcriptomic analysis results : 1. if I should only focus on uniqueMapping reads (mapped to one loci) and neglect MultiMapping read (mapped to several loci or variants)! 2. or take them all together
So my questions are :
1- Considering unique and multi-aligned reads, wouldn't this be a source of bias in my results interpretation?
2- By hiding the multi-mapped reads, could I lose information on the expression of important genes?
thank you in advance

There are several methods/tools that take into account multi-mappers and deal with them in reasonable ways (e.g. pseudo-alignment methods like salmon and more traditional alignment-based methods like RSEM). If you completely ignore multi-mappers, then yes, you will be losing information, which may or may not be valuable to you.
Both of those methods utilize the expectation-maximization algorithm (or a variant thereof) to determine the relative contributions of multi-mapping reads to the various mapping locations. With salmon, you can forego alignment completely and end up with very accurate gene quantifications in hand in like 20 minutes if you want.


How do I change the identity of a sample from spatial 'Visium' data preprocessing in Seurat v3? I'm wanting to create a merged object in Seurat using 2 10x Visium 'slices'. However, when I create an object, Seurat assigns an identity "SeuratProject" to the objects (by default I'm guessing; see attached plot). This will be problematic when I attempt to merge the objects as there will be unintended duplicates of the objects.
Here's what I'm attempting to use to create the 10x Spatial Object:
brain1 <-Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)
I'd need the ability to assign unique ident to each of these (and any subsequent) slices/data.
Thank you!!

You can re-assign orig.ident after creating the object. For example:
brain1 <-Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)
brain1$orig.ident <- "Visium"


What type of Protein ID is this? I got a FASTA file from NCBI RefSeq.
These proteins are annotated with IDs like this : protein_id=XP_019970915.1.
What type of ID is this? I would like to have them converted in TrEMBL IDs in order to do, later, a GO term enrichment but I have to know where I am starting from.
Thank you in advance

I'm not sure I understand the question, but XP_019970915.1 is an 'NCBI Reference Sequence' ID for 'erythrocyte binding antigen-181' (https://www.ncbi.nlm.nih.gov/protein/XP_019970915.1/). If you look that protein up in UniProt you can find other identifiers/info and GO annotations etc (https://www.uniprot.org/uniprot/A0A151LVZ9). Does that answer your question?


Why does the SARS-Cov2 coronavirus genome end in aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa (33 a's)? The SARS-Cov2 coronavirus's genome was released, and is now available on Genbank.  Looking at it...

    1 attaaaggtt tataccttcc caggtaacaa accaaccaac tttcgatctc ttgtagatct
   61 gttctctaaa cgaactttaa aatctgtgtg gctgtcactc ggctgcatgc ttagtgcact
  121 cacgcagtat aattaataac taattactgt cgttgacagg acacgagtaa ctcgtctatc
  ...
29761 acagtgaaca atgctaggga gagctgccta tatggaagag ccctaatgtg taaaattaat
29821 tttagtagtg ctatccccat gtgattttaa tagcttctta ggagaatgac aaaaaaaaaa
29881 aaaaaaaaaa aaaaaaaaaa aaa

  Wuhan seafood market pneumonia virus isolate Wuhan-Hu-1, complete genome, Genbank

Geeze, that's a lot of a nucleotides---I don't think that's just random.  I would guess that it's either an artifact of the sequencing process, or there is some underlying biological reason.

Question: Why does the SARS-Cov2 coronavirus genome end in 33 a's?

Good observation! The 3' poly(A) tail is actually a very common feature of positive-strand RNA viruses, including coronaviruses and picornaviruses.
For coronaviruses in particular, we know that the poly(A) tail is required for replication, functioning in conjunction with the 3' untranslated region (UTR) as a cis-acting signal for negative strand synthesis and attachment to the ribosome during translation. Mutants lacking the poly(A) tail are severely compromised in replication. Jeannie Spagnolo and Brenda Hogue report:

The 3′ poly (A) tail plays an important, but as yet undefined role in Coronavirus genome replication. To further examine the requirement for the Coronavirus poly(A) tail, we created truncated poly(A) mutant defective interfering (DI) RNAs and observed the effects on replication. Bovine Coronavirus (BCV) and mouse hepatitis Coronavirus A59 (MHV-A59) DI RNAs with tails of 5 or 10 A residues were replicated, albeit at delayed kinetics as compared to DI RNAs with wild type tail lengths (>50 A residues). A BCV DI RNA lacking a poly(A) tail was unable to replicate; however, a MHV DI lacking a tail did replicate following multiple virus passages. Poly(A) tail extension/repair was concurrent with robust replication of the tail mutants. Binding of the host factor poly(A)- binding protein (PABP) appeared to correlate with the ability of DI RNAs to be replicated. Poly(A) tail mutants that were compromised for replication, or that were unable to replicate at all exhibited less in vitro PABP interaction. The data support the importance of the poly(A) tail in Coronavirus replication and further delineate the minimal requirements for viral genome propagation.
Spagnolo J.F., Hogue B.G. (2001) Requirement of the Poly(A) Tail in Coronavirus Genome Replication. In: Lavi E., Weiss S.R., Hingley S.T. (eds) The Nidoviruses. Advances in Experimental Medicine and Biology, vol 494. Springer, Boston, MA

Yu-Hui Peng et al. also report that the length of the poly(A) tail is regulated during infection:

Similar to eukaryotic mRNA, the positive-strand coronavirus genome of ~30 kilobases is 5’-capped and 3’-polyadenylated. It has been demonstrated that the length of the coronaviral poly(A) tail is not static but regulated during infection; however, little is known regarding the factors involved in coronaviral polyadenylation and its regulation. Here, we show that during infection, the level of coronavirus poly(A) tail lengthening depends on the initial length upon infection and that the minimum length to initiate lengthening may lie between 5 and 9 nucleotides. By mutagenesis analysis, it was found that (i) the hexamer AGUAAA and poly(A) tail are two important elements responsible for synthesis of the coronavirus poly(A) tail and may function in concert to accomplish polyadenylation and (ii) the function of the hexamer AGUAAA in coronaviral polyadenylation is position dependent. Based on these findings, we propose a process for how the coronaviral poly(A) tail is synthesized and undergoes variation. Our results provide the first genetic evidence to gain insight into coronaviral polyadenylation.
Peng Y-H, Lin C-H, Lin C-N, Lo C-Y, Tsai T-L, Wu H-Y (2016) Characterization of the Role of Hexamer AGUAAA and Poly(A) Tail in Coronavirus Polyadenylation. PLoS ONE 11(10): e0165077

This builds upon prior work by Hung-Yi Wu et al, which showed that the coronaviral 3' poly(A) tail is approximately 65 nucleotides in length in both genomic and sgmRNAs at peak viral RNA synthesis, and also observed that the precise length varied throughout infection. Most interestingly, they report:

Functional analyses of poly(A) tail length on specific viral RNA species, furthermore, revealed that translation, in vivo, of RNAs with the longer poly(A) tail was enhanced over those with the shorter poly(A). Although the mechanisms by which the tail lengths vary is unknown, experimental results together suggest that the length of the poly(A) and poly(U) tails is regulated. One potential function of regulated poly(A) tail length might be that for the coronavirus genome a longer poly(A) favors translation. The regulation of coronavirus translation by poly(A) tail length resembles that during embryonal development suggesting there may be mechanistic parallels.
Wu HY, Ke TY, Liao WY, Chang NY. Regulation of coronaviral poly(A) tail length during infection. PLoS One. 2013;8(7):e70548. Published 2013 Jul 29. doi:10.1371/journal.pone.0070548

It's also worth pointing out that poly(A) tails at the 3' end of RNA are not an unusual feature of viruses. Eukaryotic mRNA almost always contains poly(A) tails, which are added post-transcriptionally in a process known as polyadenylation. It should not therefore be surprising that positive-strand RNA viruses would have poly(A) tails as well. In eukaryotic mRNA, the central sequence motif for identifying a polyadenylation region is AAUAAA, identified way back in the 1970s, with more recent research confirming its ubiquity. Proudfoot 2011 is a nice review article on poly(A) signals in eukaryotic mRNA.


Overlap in Paired-end Reads for Sequencing? Regarding Sequencing:

Do/can the paired-end reads have overlaps?

Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.


Can ortholog genes have different names among different species? To be specific, we have different lineages of E. coli. I wondered if the ortholog genes among them should all have the same name in all lineages or can they also have different names?

Gene naming is pretty subjective and even a single gene can have many names. If you're lucky orthologs will be named the same in different lineages, but you can't rely on it. Better to use gene ids (e.g. from Ensembl) and an ortholog db to map the orthologs. However, if you're looking at different strains, this is more tricky since most db don't go below the species level


Comparison of gene set enrichment statistics I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.

As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.

I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply.

I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:

library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")

# Create a design matrix and
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)),
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
##
## $correlation
## [1] -0.02036983
##
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
##
## $correlation
## [1] 0.3024643
##
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
##
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance.

As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:

               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

  How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value?

I am unsure how to interpret that comparison. It would be better to explore this  with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications.

  Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used inter.gene.cor = NA in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower.

Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (cameraPR) because it doesn't have the same results as in camera and it doesn't allow to estimate the intercorrelation.


using python to write bioinformatics pipelines tutorial I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example

use a aligner (say hisat)
get the output and process it using samtools

I was able to use subprocess from python2.7 for this purpose using samtools but i am not able to link both the processes.i.e given path(which I can use argparse) for directory with fastq files output would be processed bam.

sample code for samtools sam to bam :

import subprocess
subprocess.run(['samtools','view', '-bS',
                '../some.sam',
                '>',
                '../some.bam'])

Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and  runs it robustly.

Snakemake is my favourite. It's (nearly) pure Python and can generate reports.
Nextflow is growing in popularity and is pretty straightforward
Ruffus used to be reasonably popular, seemed fine when I used it
Bpipe is for bioinformatics
Airflow is a more industrial "big" solution but in wide use
See this big list of pipeline systems

I'm sure you can find something there.


How to initiate RNA-Seq analysis of TCGA files? I want to compare RNA-Seq datasets obtained from the TCGA to investigate how my gene of interest is implicated in different types of cancer.

I'm trying to download the data from the GDC Data Portal (https://portal.gdc.cancer.gov), but it is not user-friendly.

Can somebody give me some initial directions? I found a bioconductor package that facilitates access to the TCGA (https://bioconductor.org/packages/release/bioc/html/TCGAbiolinks.html), but first of all, I would like to understand how the files are organized, the nomenclature and how to inspect them.

thank you very much, Fabiano

The GDC API user guide (https://docs.gdc.cancer.gov/API/PDF/API_UG.pdf) provides some guidance on the structure of GDC data. So for getting RNA-Seq data from GDC, you would provide this filter to get the 551 text files with workflow_type set as "HTSeq - FPKM-UQ" and "TCGA-LUSC" for the project. There are other workflow types as well:

cases.project.program.name in ["TCGA"]
and cases.project.project_id in ["TCGA-LUSC"]
and files.access in ["open"]
and files.analysis.workflow_type in ["HTSeq - FPKM-UQ"]
and files.data_format in ["txt"]
and files.data_type in ["Gene Expression Quantification"]
and files.experimental_strategy in ["RNA-Seq"]

(You don’t need all those filters to get the 551 files; I included them to show what is available.) The txt files are returned as gzipped tar files that hold text files with rows of the form:

ENSG00000158486.12  5341.21601153

The above filters provide all the FPKM-UQ files for TCGA-LUSC. If you wanted to just get the data for normal samples as controls, that would require filtering on sample type as well, adding in the additional filter:

and cases.samples.sample_type in [“solid tissue normal”]

That will return 49 files. If you instead wished to get only the tumor files, the additional filter would instead be:

and cases.samples.sample_type in [“primary tumor”]

and that returns 502 files. Note that there are several different sample types for tumor and normal tissues; those just happen to be the ones used here.
(Note that if you are using the GDC Repository portal, adding that additional filter requires you to select “Add a Case/Biospecimen Filter” on the “Cases” filter tab.)

Since GDC is part of the NCI Cancer Research Data Commons (CRDC), instead of downloading the data, you can actually work with it in the cloud using the various CRDC Cloud Resources such as the ISB-CGC: https://www.isb-cgc.org/. In fact, all the RNA-Seq data for TCGA is in an ISB-CGC Google BigQuery table: (https://console.cloud.google.com/bigquery?p=isb-cgc&d=TCGA_hg38_data_v0&t=RNAseq_Gene_Expression&page=table)
And the above count data can be pulled using the query:

SELECT HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC"

If you want to discriminate between normal and tumor expression levels, the ISB-CGC BQ tables make it pretty easy to tag the expression values with the type of tissue involved. The sample type can be assigned by joining the previous query with the biospecimen table:

WITH bio as (
SELECT sample_barcode, sample_type_name
FROM `isb-cgc.TCGA_bioclin_v0.Biospecimen`
),
expr as (SELECT sample_barcode, HTSeq__FPKM_UQ
FROM `isb-cgc.TCGA_hg38_data_v0.RNAseq_Gene_Expression`
WHERE Ensembl_gene_id_v = "ENSG00000158486.12" AND project_short_name = "TCGA-LUSC")
SELECT bio.sample_type_name, expr.HTSeq__FPKM_UQ
FROM expr
JOIN bio
ON bio.sample_barcode = expr.sample_barcode

The GDC has listed some of the available sample types on this page: https://gdc.cancer.gov/resources-tcga-users/tcga-code-tables/sample-type-codes
Also, BigQuery has a user interface that makes exploring data sets and writing queries easy, but you need to log in with a Google ID, so the ISB-CGC has a BigQuery Table Search (https://isb-cgc.appspot.com/bq_meta_search/) that allows you to explore the available tables without logging in. If you want to run queries, Google offers a free tier that you can use to try things out and run small jobs (https://cloud.google.com/free/) and ISB-CGC offers free credits too. ISB-CGC also has examples of working in R with BigQuery within the Community Notebook Repository (https://isb-cancer-genomics-cloud.readthedocs.io/en/latest/sections/HowTos.html). A few useful notebooks would be 'How do I create cohorts of patients?', 'How do I plot a heatmap using data in BigQuery?', and 'How do I get started fast?'.


smallest length of reads in 3'quantseq? I am studying about RNA seq, and especially about 3'Quantseq(tagseq, 3primeseq). I wonder if there is a cutoff for the reads length. By this I mean that given that 3'quantseq targets the end of the genes and 3'UTR, then for example if reads are 75bp, then could be the case that only 3'UTR is mapped and not coding region? I am interested in Human Sapiens. Thank you in advance.

Simple answer: Yes, that can very well happen depending on the average fragment size, the read length and the length of the UTR.


Determining Read Groups Which Read Groups are correct:

java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84

or

java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84_S23_L005 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84_S23_L005

Or is a better way to determine read groups?

Thank you in advance,

Usually it doesn't actually matter, but if you want it to be very correct:

RGID=SomeLibraryID \
RGLB=SomeLibraryID \
RGPL=illumina \
RGPU=Some_Sequencer_ID \
RGSM=SomeSampleID

That assumes you have one library per sample and aren't splitting the output BAM files by lane.


How to safely and efficiently convert subset of bam to fastq? Question

How can I extract reads from a bam file (produced by bwa-mem) to fastq given a list of reference sequences to filter out?

Potential difficulties

maintaining FR orientation of pair end reads (in bam all the sequences are reference sequences)
keeping R1 and R2 reads
keeping quality scores in the same encoding as original fastq (default illumina phred scores in my case)
bam can be (ana usually is) sorted by coordinates

Almost there solutions

Sujai's perl solution in blobology does exact opposite - getting reads from list of references (so I could just reverse the list). The disadvantage is that the script outputs an interleaved fq file; requires unique names of mates, otherwise R1/R2 information is lost.
samtools + grep them all from fastq files

create a list of read names that do not map to filtered scaffolds. (cut will extract just read names, uniq will collapse pair end read names if they are the same). Then grep read names from fastq files and remove -- separator between hits

samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | cut -f 1 | uniq > list_of_reads_to_keep

grep -A 3 -f list_of_reads_to_keep foo_R1.fq | grep -v "^--$" > foo_R1_filtered_bash.fq
grep -A 3 -f list_of_reads_to_keep foo_R2.fq | grep -v "^--$" > foo_R1_filtered_bash.fq

filter bam & picard tools

Or I could do just the filtering part and use Picard-tools (Picard.SamToFastq), but as usual I am avoiding java as much as I can. I guess

samtools view foo.bam | grep -vf list_of_scaffols_filter \
  | java -jar picard.jar SamToFastq INPUT=/dev/stdin \
  FASTQ=foo_R1_filtered_bash.fq SECOND_END_FASTQ=foo_R2_filtered_bash.fq

The first solution does not really work for me since I do not want to rename all the reads in bam file and I want to keep the R1/R2 information (since R1 and R2 have different error profiles). Both solutions 2 and 3 I find bit clumsy and I am not sure if they are general, I might get some unexpected behaviour if one of reads is mapped second is not... They both relay on the same filtering step.

I was wondering about some pysam solution. I guess it will be much slower, but at least it will be much clearer and perhaps more general. Something like in Convert Bam File To Fasta File - there is pysam solution for fasta (not fastq), almost there...

Story

I have very fragmented reference genome. Some of scaffolds them are too small to works with, and some of them are contaminants (identified using blobtools).
I want to separate reads that are mapping to different groups to separate contaminants, short scaffolds and scaffold that will be used for downstream analysis. The reason is that if we remap all the reads to filtered reference (0.7 - 0.8 of original genome), the most of them (0.95 - 0.99) will still find a place where they map, therefore there is 0.2 - 0.3 of misplaced reads that will obviously have to affect downstream analysis, like variant calling.

This filtering idea is based logic that if the filtered duplicated genomic region will contain some small differences, they will attract their reads (and if I filter them I will improve variant calling) and if they will be exactly same, they will get reads assigned at random, so there is no harm in doing that.

I'm not aware of any pre-made program to do this, so I wrote one for you. This will take a BAM file with any ordering and produce properly ordered gzipped fastq files with the filtering as you requested. Internally, this iterates over all of the entries in the BAM file (ignoring secondary/supplemental entries and those where both mates map to your filter list), store the properly oriented sequence/quality/read name in a buffer, and then dumps that buffer entry to disk once the mate is found. This should be reasonably performant (hey, it's python, to don't expect too much), though if you happen to have indexed BAM files then one could think of ways to make this run faster.

Do check the output, since I've only run one test.


Access base aligned to particular reference position The short version: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?

The long version: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has A for SNP1 and T at SNP2, then I want to report an A,T genotype for read 1 at amplicon 1.

The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?

As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.

You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.

for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

That's not exactly correct syntax, but you get the idea.


Downloading all COI sequences from BOLD fails I have metabarcoding sequence data (COI) from bulk animal samples (including arthropoda, nematoda, annelida, mollusca) and I want to BLAST all of these sequences. I used following command to do this: blastn -remote -db nt -query COI_all.fasta -num_alignments 2 -out COI_blasted.txt. However this results in errors similar to this post: https://www.biostars.org/p/359971/ .

These errors probably appear due to the number of sequences in my file (around 700) and
the remote connection is thus interrupted.

I found that a solution would be to use blastn with a local database and since the samples are so diverse, I would like to download ALL animal COI sequences from BOLD (or gen bank). It would not be a problem if non-animal (e.g. plant) sequences would also be included.

I think the BOLD database would be great to BLAST my sequences to. However, I'm currently struggling to find a good way to download all animal COI sequences from BOLD.

When entering COI-5P as search term on http://v4.boldsystems.org/index.php/Public_SearchTerms I receive error: Your search terms resulted in too many matching terms. Please try again with more specific search criteria.. I could likely download the sequences from all the phyla etc seperately and merge them, but I'd rather just download 1 file.

I also tried to use the API by running: wget http://v4.boldsystems.org/index.php/API_Public/sequence?marker=COI-5P. A download starts but around 3.7 MB download, it is stuck and the file I receive only contains ~5000 sequences.

UPDATE: I've contacted BOLD about the stalling behavior and this is their reply: "This issue is because of the large API request that retrieves millions of records, which our system does not handle. Please break up the search by smaller groups, such as classes."

Does anyone have a solution to download all COI sequences from BOLD in one file?

I could also download COI sequences from gen bank using the ftp://ftp.ncbi.nlm.nih.gov/blast/db/ URL, but I'm not sure which exact files I need. For 16S, 18S,.. it is obvious, but not for COI. Any suggestions?

Thanks for the help.

A straightforward way of getting the sequences from GenBank is from their website using this URL:
https://www.ncbi.nlm.nih.gov/nuccore?term=COI[Gene Name]

There you can click on "Send to" -> select "Gene Features" and "FASTA Nucleotide" -> Click "Create File".

There are >3 mio hits, so the download could take a while.