Do you know any virtual screening libraries with small, soluble compounds? I had a quick question on virtual libraries. I'm hoping to perform a virtual screen using a library that contains very small, soluble compounds - think glycerol and smaller (<100MW). As far as I can tell, the majority of virtual libraries contain mostly larger, drug-like compounds. Is anyone aware of any libraries that might contain ideally 10,000+ or so of these tiny soluble compounds?
If not, any idea what the process might be (and how difficult it would be) to construct a custom virtual library enriched for these compounds?
Thank you!!

Lipinski and logP
Most datasets will be filtered by the Lipinski rule, plus in the some cases, such as Enamine Real whether they can be bought.
One criterion in the Lipinski rule filter is the logP, or octanol-water partition coefficient. This is super important, because a compound that cannot cross the membrane is likely of limited uses —as a consequence prodrugs are a commonly taken approach.
However, albeit obvious, not all datasets in Zinc are filtered by the Lipinski rule, they just need more browsing...
Natural?
However, you may have a valid reason to eschew logP.
Nature tries its best to have intermediate metabolites that have hydrophilic handles attaches, eg. CoA, phosphates and sugars. Therefore, if you want a natural metabolite this is a completely different question. Many of these are not very purchasable, but a PDB-derived dataset of natural ligands will likely be purchasable or easily extractable. Alternatively Zinc has several natural catalogues...
Too small
Glycerol is a 3 carbon molecule and has a MW of 90 Da. Many of the old drugs are small, but not that small. For example aspirin is 180 Da, which is twice the size of glycerol. Lipinski rule says not to go over 500 Da and 250-400 Da is a lovely range. Say you have a perfectly bound glycerol, with a perfect ligand efficiency of -1.5 kcal/mol/atom, that would mean it would bind at –9 kcal/mol in the best case scenario, i.e. it is too small to bind anything with nanomolar affinity and will actually bind a lot of things. Glycerol after all is commonly cocrystallised in protein...
The GDB datasets have enumerated compounds of a given size, but don't go that low.
Docking
There will be two further complications that will hinder the accuracy.

protonated and deprotonated atoms
polarisable groups

The former requires correct conformers, the latter requires a force field that has a Drude particle, e.g. in openMM.
Caveat against sorting by binding score
Also, it should be said that if you sort by binding score (kcal/mol or kJ/mol) you will get bigger compounds (hence the whole point of ligand efficiency) and compounds with high logP. The former is because they will have more contact points to make, the latter is because given that ∆∆G is ∆G_bound - ∆G_unbound, an unbound high-logP compound is highly unhappy in aqueous solution.


How to convert featureCounts to FPKM? I have seen many posts regarding counts to RPKM and TPM. I haven't seen any post for counts to FPKM.

I have RNA-Seq data which is paired-end reads. Extracted the counts using featureCounts for all the samples.

There is a function to convert counts to RPKM: using the gene_length

rpkm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(counts) * 1e6
}

I know that RPKM is mainly used for single-end reads data. Do you think I can use the above function for converting counts to FPKM as my data is paired end? [TCGA data: HTSeq counts are converted to HTSeq FPKM]

Somewhere I have seen that the same function applies for single end data which will be RPKM and also for paired end data which will be FPKM. Am Ir right or wrong?

If not can anyone show some function or code to convert counts to FPKM please.

thanq

I have seen many posts regarding counts to RPKM and TPM.

There’s your answer then: FPKM = RPKM. It’s simply a more accurate name.

Speaking of RPKM for paired-end data is discouraged because the reference to “read” in this context lends itself to ambiguity. But mathematically the quantity is the same: we are counting fragments, not individual reads (of which each fragment has two, for paired-end data).

But as mentioned in the comments, there are good reasons against using FPKMs.


How to find a homolog that has a PDB structure available I tried to find a homolog that has a PDB structure available so I can use this PDB file for comparative modeling. I ran a BLAST search but none of the search results seemed to have a known structure. Is there any way to find the results I want?

You could use HMM-based comparison, for example with HHpred. This is generally the best way to find structural templates because structure is more conserved than sequence, and comparing HMMs will turn up hits of the same structure but minimal sequence overlap.

If HHpred turns up nothing you are in the region of de novo modelling, which is a different problem. However, it has advanced in recent years to the point where it may be useful to you. For example, you could try our method DMPfold.


Mark BWA-SW split alignments in output for long reads Is there a way to remove split alignments of single-end long reads from BWA-SW output? BWA-MEM has an option to include or flag split alignments in the result, but it seems BWA-SW method include them in the result unless the split occur in paired-end reads.

From BWA manual, it is advised to use BWA-MEM in your longs read case (unless you have a good reason not to):

  BWA-MEM and BWA-SW share similar features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and more accurate.

I guess that BWA-SW returns several lines per transcript (based on local alignments). In that case, the issue has already been reported. The user asks:

  How can I differentiate between a read that has multiple hits in the reference and a read that has a split alignment?

The author replies:

  In principle, it is hard to tell. Bwa-mem makes an arbitrary choice.
  Mappings flagged with 0x100 are multiple hits. Mappings flagged with
  0x800 are from a chimeric alignment.

So you could actually filter by flag. But I am afraid the output could include something you do not necessarily want to (since there is no specific flag for the split reads).

Maybe a better solution: you could switch to some mapper that will report the proper split in the hit of the read map, such as minimap2 from the same author. It is suitable for long reads and comes with a splice -x feature.

Then, the laziest solution I see here would be to remove the split alignments exploiting the sam format and in particular the CIGAR string. Meaning, to remove the sam fields containing a N in their CIGAR (aka, skipped region from the reference). Something like this would work:

cat alignment.sam | awk '$6 !~ /N/ {print $0}'

Hope this can be somewhat helpful.


Bioinformatics flowchart for Materials and Methods I am going to start work on the Materials and Methods for a paper I contributed to a lot in the bioinformatics processing part. I have seen some papers nowadays tend to first present the processing steps that they followed as a flow chart, to show the reader what happened generally before going into details in the written M&M section.

Are there any specific guidelines for creating this flowchart? Do you have any tips or advice from personal experience?

(I am aware that this may be off-topic. I have a meta question about it here)

If you're using a workflow management system for pipelining like snakemake or ruffus, you can actually build a flowchart based on the rules/functions in your code:

Snakemake example with the --dag flag:

Ruffus example using pipeline_printout_graph:

You'll probably still need edit the image aftwards in software like inkscape or illustrator though.


Obtain CD (cluster of differentiation) marker expression information There are 394 CD marker genes. Various resources (BD handbook, Abcam guide) contain information about which cell types express which CD markers. What is the best way to obtain or download this data so that it can be used with programs such as R?

How can we download similar annotation data on expression for any gene or cell type?

I can offer this as a way to scrape the data from the Abcam resource:

library(rvest)
library(dplyr)

antigens <- read_html("http://www.abcam.com/primary-antibodies/human-cd-antigen-guide") %>%
  html_node("table") %>%
  html_table(header = TRUE)

# what targets expressed in "Langerhans" ?
antigens %>%
  filter(grepl("Langerhans", `Cellular Expression`)) %>%
  select(`Target Name`)

   Target Name
1         CD1a
2         CD1b
3         CD1c
4         CD1e
5         CD14
6         CD39
7         CD50
8         CD69
9         CD74
10        CD83
11       CD101
12       CD206
13       CD207

The data could be further processed by, for example, splitting the Cellular Expression column on either and or , (comma) to try and make one row per cell type.


Is there a way to tell which chromosome a gene is on, by looking at the "Chromosome/scaffold name" I recently got a data set, from which I need to figure out which chromosome a gene is from, but the head of the data reads like:

Gene ID Description     Gene type       Gene End (bp)   Gene Start (bp) Strand  Associated Gene Name    Chromosome/scaffold name

while one entry of the data is like:

ENSG00000252303 RNA, U6 small nuclear 280, pseudogene [Source:HGNC Symbol;Acc:HGNC:47243]       snRNA   67546754        67546651        1       RNU6-280P       CHR_HG2128_PATCH

It looks like the last column is the information about which chromosome the gene comes from, but how to interpret the "CHR_HG2128_PATCH"?

What you are looking at is a patch, not an actual chromosome location. HG-2128 is an issue ID, something that the Genome Reference Consortium uses to track issue with references. Information regarding that specific patch is here. The issue is described as a "Possible misassembly or indel variation in GRCh38 within AL049860.8."

During the alignment step before you received your data, someone used a top level assembly version of the reference genome. This includes extra, accessioned scaffolds tacked on to the end of the file that contain the fix. These scaffolds represent what will be found in the next major release of the genome. The fixes aren't applied directly to the reference as soon as they are found due to the inevitable disruption of coordinates. So reads mapped to this FIX patch rather than the primary original scaffold in the primary assembly, on which the RNU6-280P gene resides.

The coordinates are slightly different. In the current GRCh38.p12
 assembly the RNU6-280P gene is located at chr6:67546651-67546754 (source). You can see that the coordinates are very close to those indicated in the table you provided (Gene Start, Gene End) but are indeed different.

Here is some more information on patches.


How to enlarge a section of y-axis? I am using ggplot2 to plot a scatter plot.

library(ggforce) # required by facet_zoom

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

As most of my data points are in the range of 0 - 200 on y-axis, I want to zoom this section. I tried doing it using facet_zoom from ggforce package :

g1 + facet_zoom(y = Average_Gene_Expression > 0 & Average_Gene_Expression < 200)

Instead of two plots in one figure, I want one plot (like first figure), where 80% of y-axis displays the data points in the range of 0 - 200 and 20% of the y-axis displays the remaining and x-axis unaltered. Which function/package should I use do it in R?

I would strongly discourage you from making discontinuous axis, it's going to be very confusing for a reader.

The facet plot you proposed seems like a good solution to me. Alternatively you can use log transformation. To demonstrate I made it on simulated data that look appox like yours :

set.seed(940401)

data.plot <- data.frame(Methyl_Average = c(0.05, rbeta(699, 2, 5) + 0.15, runif(300, 0.18, 0.85)),
                        Average_Gene_Expression = c(165, exp(rnorm(999, 3))))

library(ggplot2)
require(gridExtra)

g <- ggplot(data.plot, aes(x = Methyl_Average, y = Average_Gene_Expression))
g1 <-  g + geom_point(color = "steelblue") + theme_bw(base_family = "Times") +
       ggtitle( "Scatter Plot") +
       theme_bw(base_family = "Times") +
       theme(plot.title = element_text(hjust = 0.5))

sc_plot <- g1 +
           labs(x = "Differential Methylation", y = "Differential Expression")
logsc_plot <- g1 +
              labs(x = "Differential Methylation", y = "Differential Expression (log10 scale)") +
              scale_y_continuous(trans='log10')

grid.arrange(sc_plot, logsc_plot, ncol = 2)


Gene information from chromosome location I have a file that contains species name, chromosome, start and end locations like the following format
hg19,chr8,77778733,77779026. I would like to know some information about the regions 77778733-77779026. How can I quickly find out if it lies in an intron or exon of a gene? Or it is in the upstream or downstream regions of a gene? Is there a tool or web service from where I can get this information?

How can I quickly find out if it lies in an intron or exon of a gene?

I have an answer on Biostars that might help: https://www.biostars.org/p/124515/#124522

  Or it is in the upstream or downstream regions of a gene?

One option is to build a table of annotated bins and then map your regions-of-interest to those annotations.

These annotations would consist of whether the bin is inside a gene, or, after picking the nearest gene by strand, locate the bin to an upstream region of the gene TSS or to an downstream region of the gene TES.

To start, one could download stranded gene annotations:

$ GENCODE_URL=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.basic.annotation.gff3.gz
$ GENCODE_GENES=${PWD}/../data/gencode.v28.genes.bed
$ wget -qO- "${GENCODE_URL}" | gunzip -c | awk '($3=="gene")' | gff2bed - | awk -v FS="\t" -v OFS="\t" '{ n=split($10,a,";"); p=a[4]; n=split(p,a,"="); $4=a[2]; print $0; }' | cut -f1-6 > ${GENCODE_GENES}

Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):

$ ASSEMBLY=hg38
$ GENOME=${PWD}/../data/${ASSEMBLY}.bed
$ CHOP_SIZE=1000
$ CHOPPED_GENOME=${PWD}/../data/chopped.${ASSEMBLY}.bed
$ fetchChromSizes ${ASSEMBLY} | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | grep -ve "\.*_\.*" | grep -v "chrM" | sort-bed - > ${GENOME}
$ bedops --chop ${CHOP_SIZE} ${GENOME} > ${CHOPPED_GENOME}

Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories.

Here's a script I wrote which uses the output of closest-features to do that. This could be easily modified for any set of distance categories:

#!/usr/bin/env python

import sys
import numpy as np

'''
Input should be bin midpoints run through closest-elements, against some set of BED6+ (stranded) gene annotations
'''

def assign_to_category(dist):
  if dist < -50000:
    return '-50kb_to_-Inf'
  elif dist < -20000 and dist >= -50000:
    return '-20kb_to_-50kb'
  elif dist < -10000 and dist >= -20000:
    return '-10kb_to_-20kb'
  elif dist < -5000 and dist >= -10000:
    return '-5kb_to_-10kb'
  elif dist < -1000 and dist >= -5000:
    return '-1kb_to_-5kb'
  elif dist < 0 and dist >= -1000:
    return 'TSS_to_-1kb'
  elif dist == 0:
    return 'Gene'
  elif dist > 0 and dist <= 1000:
    return 'TES_to_+1kb'
  elif dist > 1000 and dist <= 5000:
    return '+1kb_to_+5kb'
  elif dist > 5000 and dist <= 10000:
    return '+5kb_to_+10kb'
  elif dist > 10000 and dist <= 20000:
    return '+10kb_to_+20kb'
  elif dist > 20000 and dist <= 50000:
    return '+20kb_to_+50kb'
  elif dist > 50000:
    return '+50kb_to_+Inf'

for in_line in sys.stdin:
  (bin_str, upstream_gene_str, upstream_gene_dist, downstream_gene_str, downstream_gene_dist) = in_line.rstrip().split('|')
  try:
    upstream_gene_dist = int(upstream_gene_dist)
  except ValueError:
    upstream_gene_dist = np.nan
  try:
    downstream_gene_dist = int(downstream_gene_dist)
  except ValueError:
    downstream_gene_dist = np.nan
  upstream_gene_elems = upstream_gene_str.split('\t')
  downstream_gene_elems = downstream_gene_str.split('\t')

  '''
  We have six possible scenarios, depending on the minimum distance and strandedness of closest elements
  '''

  '''
  1) If either distance value is NaN, the bin is assigned to the non-NaN distance element
  '''
  if upstream_gene_dist is np.nan and downstream_gene_dist is np.nan:
    raise SystemError("Bad NaN input on:\n[{}]\n".format(in_line.rstrip()))
  elif upstream_gene_dist is np.nan:
    gene_id = downstream_gene_elems[3]
    if downstream_gene_elems[5] == '+':
      # bin is closer to downstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      sys.stdout.write('{}\t{}\t\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue
  elif downstream_gene_dist is np.nan:
    gene_id = upstream_gene_elems[3]
    if upstream_gene_elems[5] == '+':
      # bin is closer to upstream TES
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to upstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    continue

  '''
  2) If either distance value is zero, the bin is assigned to "Gene" to state overlap with a gene
  '''
  if upstream_gene_dist == 0:
    gene_id = upstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  if downstream_gene_dist == 0:
    gene_id = downstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  '''
  3) If the upstream strand is positive, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  4) If the upstream strand is negative, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TSS
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  5) If the upstream strand is positive, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  6) If the upstream strand is negative, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  Sanity check: If we get here, something went wrong
  '''
  raise SystemError("Bad logic on:\n[{}]\n".format(in_line.rstrip()))

To use it, e.g.:

$ MP_LEFT=500
$ MP_RIGHT=-499
$ MP_LEFT_UNDO=-500
$ MP_RIGHT_UNDO=499
$ CATEGORIZED_BINS=${PWD}/../data/categorized_bins.${ASSEMBLY}.bed
$ ./assign_chops.py < <(bedops --range ${MP_LEFT}:${MP_RIGHT} --everything ${CHOPPED_GENOME} | closest-features --dist - ${GENCODE_GENES}) | sort-bed - | bedops --range ${MP_LEFT_UNDO}:${MP_RIGHT_UNDO} --everything - > ${CATEGORIZED_BINS}

The output from this file annotates bins with the closest feature to a gene. It looks something like this:

$ more categorized_bins.hg38.bed
...
chr1  32000 33000 TES_to_+1kb RP11-34P13.3
chr1  33000 34000 TES_to_+1kb FAM138A
chr1  34000 35000 Gene  FAM138A
chr1  35000 36000 Gene  FAM138A
chr1  36000 37000 Gene  FAM138A
chr1  37000 38000 TSS_to_-1kb FAM138A
chr1  38000 39000 -1kb_to_-5kb  FAM138A
...

For example, the region chr1:33000-34000 is closest to the 1kb window downstream of the TES of FAM138A.

Once you have this file, you can quickly map it to your sorted regions-of-interest file:

$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed > answer.bed

The file answer.bed will contain your regions of interest and the annotated bins which overlap them.


Tool for predicting interactions in the cell What tools are available to predict, based on the structure of a certain protein, its interactions within a cell?

For example, I am considering the split GFP protein and I am trying to predict if there will be any interactions which might hinder its binding efficiency to the other half of the GFP protein.

Note that I am not only looking for protein interactions, but other molecules as well.

You could try one of these tools to predict protein-protein interactions:

Struct2Net

Given two protein sequences, the structure-based interaction prediction technique threads these two sequences to all the protein complexes in the PDB and then chooses the best potential match. Based on this match, the method generates alignment scores, z-scores, and an interfacial energy for the sequence pair. Logistic regression is then used to evaluate whether a set of scores corresponds to an interaction or not. The algorithm is also extended to find all potential partners given a single protein sequence. Further details about the method are described here.

PRED

Pred_PPI is a web-based system that serves for predicting PPIs from different organisms. This server is freely available to any researcher wishing to use it for non-commercial purposes. Based on auto covariance (AC) and support vector machine (SVM), this tool is capable of predicting PPIs for any target protein pair only using their primary sequences, and assigning an interaction probability to each SVM prediction as well. So the user can use this tool to predict novel PPIs with high confidence.

For protein-RNA interactions, you can try:

catRapid

Through the calculation of secondary structure, hydrogen bonding and van der Waals contributions, catRAPID is able predict protein-RNA interaction propensities with great accuracy (up to 89% on the ncRNA-protein interaction database, NPinter).


RNAseq data analysis I am currently doing a Differential Gene Expression of RNAseq data in Lung Adenocarcinoma using TCGAbiolinks. In the data Preprocessing step TCGAanalyze_Normalisation, I am confused as to which method to select for normalization (geneLength or gcContent). What are the possible implications of both methods and how to decide which is the best suited model for the analysis?

Is there a better way to perform the DGE ?

The TCGAbiolinks reference manual provides a clear explanation of the two methods:

  Normalization for RNA-Seq Numerical and graphical summaries of RNA-Seq
  read data. Withinlane normalization procedures to adjust for
  GC-content effect (or other gene-level effects) on read counts: loess
  robust local regression, global-scaling, and full-quantile
  normalization (Risso et al., 2011). Between-lane normalization
  procedures to adjust for distributional differences between lanes
  (e.g., sequencing depth): global-scaling and full-quantile
  normalization (Bullard et al., 2010).

Ideally, you would know if your samples spanned different lanes or not, which would indicate which you should use - "geneLength" if so, "gcContent" if not. If you're using TCGA data, that metadata is likely available somewhere.


Autodetect max number of cores and pass as an argument in Nextflow I am creating a pipeline in Nextflow. One step is creating a pangenome with Roary. Roary takes threads as an argument and if no number of threads is supplied as an argument it defaults to one.

Is there a way in Nextflow to pass the maximum number of available threads as an argument? At the moment the best I can come up with is defining an arbitrarily large number of threads as an argument which will cap out and run with the max number of threads. I would like to know, is there a 'cleaner' or more canonical way of doing this?

At the moment my process looks like this:

process roary {

    publishDir "${params.outdir}/roary", mode: 'copy'

    cpus 200

    input:
    file gff from gff.collect()

    output:
    file("*") into roary
    file("pan_genome_reference.fa") into pan_genome
    file("pan_genome_sequences/*") into alignment_files
    file("gene_presence_absence.Rtab") into gene_presence
    set file("*accessory*"),
        file("*.Rtab"),
        file("_*"),
        file("*.txt"),
        file("*csv") into roary_out

    script:
    """
    roary -p ${task.cpus} -e -n -v -z $gff
    """
}

Thanks in advance. Also, apologies as I can see why this might not strictly be bioinformatics, but it's where Nextflow has the most traction.

I have not read the full Nextflow manual in depth, but I think there is no Nextflow command to detect the maximum number of threads you can use.

However, you can try to automatize this task by using the nproc command.

nproc gives you the number of processing units available. An idea would be to include something such nproc - 1 as the max number of threads for Roary.


How can I specify to DEseq2 to only perform comparisons on which I am interested with? I am currently performing a large RNA-seq analysis from mice PBMCs. The dataset contains around 6,000 transcriptomic profiles and I would like to use DESeq2 to identify the sets of differentially expressed genes in the different conditions. In total, I have 100 biological stimulations, and for each stimulation I have 30 control samples and 30 samples treated with a molecule of interests (so is have: (30 controls+30 treated) * 100 stimulations = 6,000 samples).

I want to identify, for each stimulation, the sets of differentially expressed genes between control samples and treated samples. I do not want to compare samples from the different stimulations. In total, I would like to have thus 100 lists of differentially expressed genes.

I have started to use deseq2 to identify these lists but deseq2 is spending lot of time to perform comparisons on which I am not interested with (comparisons between the biological stimulations).

For now, I have a table sampleTable which looks like that:

I am using DEseq2 using the following command:

DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)

Could you please help me with that? How can I specify to deseq2 to not perform the comparisons between the biological stimulations, but rather between control and treated sample within each stimulation ?

Thank you and best,

You can specify the exact comparisons you want in the results() function. So:

dds = DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds = DESeq(dds)
res = results(dds, contrast("condition", "treatment1", "control1"))

Note that condition should not be, "stim001-control1, stim001-control2, etc.", but instead, "stim001control, stim001control, stim001control, stim001treatment, etc.". Don't put minus signs in your levels and don't number your samples in them (they're not in different groups.

The last command would be repeated for each of the comparisons. Note that this will still be quite slow due to the number of samples and the number of groups. With so many samples, you might just make a separate sampleTable for each of the comparisons you want to make. Each would then only contain the samples relevant for that comparison, so 60 total in each.

Alternatively, if you want to fit all of the groups at once then it's likely that limma/voom will prove to have better performance. It uses different math that happens to be quicker with large models.

BTW, if the slowness is coming from DESeqDataSetFromHTSeqCount() then you can mere everything on the command line into a single matrix. A couple lines of python or even the join command can do that quickly enough.


How to plot circular chromosomes with Gviz? I am using Gviz to plot several tracks. I create a GenomeAxisTrack() and then several AnnotationTrack()s that I feed to plotTracks(). This produces a stacked output of tracks with annotated arrows as I want but the displays are linear. In some cases I have the complete genomes and I know they are circular (like some bacterial and viral chromosomes), how can I display the plots as circular?

I tried feeding to AnnotationTrack(); a GRanges() object with a Seqinfo() object specifying the isCircular=TRUE flag, but this didn't change the resulting plot.

Here is minimal code to create a linear plot, I would like one where the first and last arrow join up in a circle:

library(Gviz)
library(stringr)

track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('+','+','-'), stacking="dense", fill=c('red','blue', 'black'))
png(file=paste("test","png", sep="."))
plotTracks(track)
dev.off()

Note, I am not interested in alternative software recommendations, I would like to solve this with gviz if it is possible. Also asked this here.

The primary author of Gviz confirmed to me via email that this is not possible:

  Gviz will only give you the traditional linear tracks. There are a
  bunch of solutions available for visualizing circular genomes, most of
  which are more interactive in nature and not implemented in the rather
  static R graphics context. GGplot2 has some support for circular axes,
  too. If you want to stay within the R environment I’d start there.
  Ggbio is another option providing a somewhat higher-level interface.


python package for NNI neighbors I am working on protein sequence data files for reconstructing phylogenetic tree and I need to generate all NNI-neighbours of the tree (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation). Is there any python package/sorce code I can access to, which can retrieve all NNI neighbors of a tree.
Earlier I have used biopython 1.72 to generate tree from my dataset. But for now I'm looking for a python package which can do NNI search on my tree and return all NNI-neighbors.

Both the ETE Toolkit and DendroPy offer methods of node removal and reattachment - I would think it would be easy enough to hack together an NNI movement function with them.

You could also investigate the source code of BioPython's NNITreeSearcher() function.


storing SNPs in the genotype file for each chromosome in a separate file I have already asked this question in another forum but have not got an appropriate answer, so wondering if anyone here can help me?
I have a big file with over 3 million columns (SNPs genotypes) and 800 rows (one row for each individual).

The first row is the SNPs_ID name with the SNP position on each chromosome (see the example below please), in total there are 17 chromosomes. I want to extract SNPs for each chromosome and store them in a separate file.
I have done it with "R", but I am looking for a faster solution and I am wondering if anyone can help me to figure out another solution with bash, perl or python codes?

here is an example of my input file (SNPs on chromosome 1,2,3,4...,17)

"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"....... "Ha412HOChr17:303402"
0 1 0 0 ......0
0 1 0 0 ......0
0 1 0 0 ......0
0 2 0 0 ......0
0 1 1 1 ......1
0 2 0 0 ......0

and this is an example of my desired output for example for chromosome 1:

out.chrom1
"Ha412HOChr01:180159" "Ha412HOChr01:210724" "Ha412HOChr01:303270" "Ha412HOChr01:303280"
0 1 0 0
0 1 0 0
0 1 0 0
0 2 0 0
0 1 1 1
0 2 0 0
0 0 0 0
0 2 0 0
0 1 2 2

I believe a pandas approach would be fast enough (probably faster than any other Python approach without pandas).

After converting your table to a pandas data frame, you can do something like df.filter(regex=(".*Chr01")) for the first chromosome then write the resulting data to file with the to_csv() method.

import pandas as pd

data = pd.read_csv('input.txt', delimiter = "\t", index_col = False)

for i in ["Chr01", "Chr02", "Chr03"]:
    regex_pattern = ".*" + i
    print(data.filter(regex=regex_pattern))


Axes for ggplots can be controlled via xlim() and ylim() functions. You can try:

your_ggplot_object +
  ylim(-5,20) # assuming that you want the gene labels to be readabe, not clear from your post


Telling grep to treat N as [ATCG] Okay so I'm using grep to try and get a preview of some trimming operations that are not going as expected..
Lets say that my sequence in the FastQ file is: ATNGCNATCG
What I want to do is..
grep "ATCGCTATCG" my.fastq
..and match the sequence given above
Surely there is some way or some existing tool that I can use besides doing:
grep "[A|N][T|N][C|N][G|N][C|N]...etc."

If you want to stick to grep, use a scripting language such as Perl to generate the regex programmatically. For example:
perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' ATCGCTATCG

Prints:
[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]

You can use it in grep like so:
grep '[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]' <<< ATNGCNATCG

Prints:
ATNGCNATCG

If that works for you, you could make it into a little bash function that also runs the grep. Add these lines to your ~/.bashrc:
grepN(){
    seq="$1"
    file="$2"
    pattern=$(perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' "$seq")
    grep "$pattern" "$file"
}

You can now run
grepN ATCGCTATCG my.fastq

Of course, this is not a good idea since the sequence might be in different lines, but that's what you were doing originally.


add task to nextflow and resume in test.nf (dsl2) I have
   ...
   workflow{
       a = Channel...
       processA(a)
       processB(processA.out.emitA)
       processC(processB.out.emitB)
    }

The pipeline completed successfully.
Now I want to add processC_2 that takes the output of processB.
How do I run the workflow to bypass processA execution?
I added processC_2 to test.nf and then did:
nextflow run test.nf -resume run-name

But it doesn't run the newly add process.
How do I achieve my goal?
Thank you!

It is actually very hard to answer this using only the information you have provided. It'd be necessary to actually see your test.nf file and read the processes. The execution seems right and I'd guess that the processB may be not emiting any output as shown in https://www.nextflow.io/docs/latest/process.html#outputs.
It is also possible that when using the -resume together with a specific run-name you may be recovering the state of that specific run as shown here https://www.nextflow.io/blog/2019/demystifying-nextflow-resume.html. You can try to execute it without the run-name or without the -resume flag so you can check whether the problem is with the emition and channelling of outputs between processes, or with the -resume flag.


dividing genome into non-overlapping windows using R I have nearly 10 million SNPs located on 10 chromosomes. I want to divide the genome into non-overlapping windows of 15, 20 and 30 kb. Here is part of my SNP table:

head (sap_ids)
       snp_id    chr  pos
    Chr01__15043   1 15043
    Chr01__15079   1 15079
    Chr01__15139   1 15139
    Chr01__15165   1 15165
    ...
    ...
    ...
    Chr17__214708424  17 214708424
    Chr17__214708451  17 214708451
    Chr17__214708484  17 214708484
    Chr17__214708508  17 214708508
    Chr17__214708574  17 214708574

I have been using this code, but it gives me wrong output because each chromosome starts from position 1

win_size<-c(15000,30000,50000,100000)
res<- cbind(snp_ids,
       data.frame(lapply(setNames(win_size, paste("window",win_size, sep = "_")),
       function(x)as.numeric(ceiling(snp_ids$pos/x)))))

For instance, as I you see in the example below in window 4 I get SNPs from chromosomes 1, 3 and 5:

snp_id     chrom poistion   window
Chr01__58332   1 58332            4
Chr01__58335   1 58335            4
Chr01__58341   1 58341            4
Chr01__58450   1 58450            4
Chr01__58471   1 58471            4
Chr01__58530   1 58530            4
Chr01__58542   1 58542            4
Chr01__58641   1 58641            4
Chr03__45457   3 45457            4
Chr03__45604   3 45604            4
Chr04__56873   4 56873            4
Chr04__57387   4 57387            4
Chr04__57399   4 57399            4
Chr04__57528   4 57528            4
Chr04__58419   4 58419            4
Chr04__59670   4 59670            4
Chr04__59704   4 59704            4

I guess the solution would be to run the code above by looping through each chromosome I tried something like this:

for (i in 1:10){
       res <- cbind(cats,
             data.frame(
               lapply(setNames(win_size, paste("window", win_size, sep = "_")), function(w)
                 paste(cats$chr[cats$chr == i], ceiling(cats$pos/w), sep = "_"))
               ))
       }

But still not working appropriately (column 4 should be 1_1, for the chr 1)!

Chr01__912    1  912          2_1
Chr01__944    1  944          2_1
Chr01__1107   1  1107         2_1
Chr01__1118   1  1118         2_1

Can anyone help me to figure out how can I modify the code to create these non-overlapping windows properly?

The question is a bit confusing to me, at it's core I understand you want non-overlapping windows across chromosomes.

One way to achieve is to use GenomicRanges::tileGenome function, which needs the chromosome lengths as input, e.g.

library(Biostrings)
mygenome <- readDNAStringSet(list.files(mypath,"mygenome.fa$",full=TRUE))
chrSizes <- width(mygenome)
names(chrSizes) <- names(mygenome)
print(chrSizes)
#     Chr1     Chr2     Chr3     Chr4     Chr5
# 18585056 19698289 23459830 26975502 30427671

This can than be fed into tileGenomes:

library(GenomicRanges)
bins   <- tileGenome(chrSizes, tilewidth=5e5, cut.last.tile.in.chrom=T)

This creates a GRange object which gives you a whole repertoire of post processing options, i.e. overlapping with SNPs or NGS reads (I use the latter for genome wide coverage plots):

print(bins)
# GRanges object with 240 ranges and 0 metadata columns:
#         seqnames               ranges strand
#            <Rle>            <IRanges>  <Rle>
#     [1]     Chr1   [      1,  500000]      *
#     [2]     Chr1   [ 500001, 1000000]      *
#     [3]     Chr1   [1000001, 1500000]      *
#     [4]     Chr1   [1500001, 2000000]      *
#     [5]     Chr1   [2000001, 2500000]      *
#     ...      ...                  ...    ...
#   [236]     Chr5 [28000001, 28500000]      *
#   [237]     Chr5 [28500001, 29000000]      *
#   [238]     Chr5 [29000001, 29500000]      *
#   [239]     Chr5 [29500001, 30000000]      *
#   [240]     Chr5 [30000001, 30427671]      *
#   -------
#   seqinfo: 5 sequences from an unspecified genome


Multiple sequence alignment with R I want to incorporate Multiple sequence alignment into our Shiny app, i.e. allow the possibility of selecting different nucleotide sequences, and do a multiple sequence alignment in real-time. The user would define as input different nt. sequences and then a tree will pop up in the app. Do you think this is feasible?

Any recommendations and feedback are really appreciated :)

The ‹msa› package on Bioconductor does exactly that.

It doesn’t hand the results on a silver platter, though: you’ll need to read the vignette carefully to learn how to use it. But after that it’s pretty powerful.


Extract reads from bam files by their @RG How could I extract reads from bam files by their read groups @RG ?, I've got one file containing all reads for 5 samples, the SM is appears NONE, So, I want to extract each sample reads by the @RG tags.

Typically I use samtools for operations like this.  Specifically I use samtools view with either -r or -R flag depending on the use case.

-r STR
Output alignments in read group STR [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.

-R FILE
Output alignments in read groups listed in FILE [null]. Note that records with no RG tag will also be output when using this option. This behaviour may change in a future release.


This does not look a great confusion matrix, but it is a very, very cool approach. The cool thing about machine learning is the ability to examine lots of models within confusion matrices.

To answer your question, you want a heat map of the regression weights to establish who is who, i.e. which patient is false positive etc... This is the only way I'm aware of to analyse a neural network. This looks like R, I only know Scikit learn.


plot correlation between several bam files I have several mab files from medip-seq. they supposed to be replicates and I have to draw their correlation. First I used multiBamSummary from deeptools to merge all the bam files based on bi/bed. Now when I am using the plotCorrelation I do not get any file. the argument is :

plotCorrelation --corData myMultiBamSumary.npz --corMethod spearman --whatToPlot scatterplot --skipZeros --removeOutliers --plotFile Myplot --plotFileFormat pdf

Is that the correct strategy to pot the correlation? Is there any suggestion?

Thanks alot

That's a completely normal way to plot correlation between them. Spearman correlation was a wise choice on your part, since it's rank based and you don't have to care as much about normalization.

Update: It turns out that matplotlib is broken and I need to have a new deepTools release to work around that. That should be available sometime Sunday.


How to get results from Homo.sapiens package in bioconductor for a specific reference I want to use the Homo.sapiens package in Bioconductor to retrieve the chromosome location start and end for each gene symbol in a specific reference (e.g. hg19 or hg38).

Right now, I am using the following code:

> select(Homo.sapiens, keys=c('KRAS'), cols <-c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype='SYMBOL')

  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863

My question is two-fold:

How can I specify the human genome reference I want to retrieve the data from? (e.g. hg19 or hg38)
Why do I receive multiple and different mappings through Homo.sapiens for some genes? If I check the locations through GeneCards I get:

chr12:25,204,789-25,250,936 (GRCh38/hg38)
chr12:25,357,723-25,403,870 (GRCh37/hg19)

I couldn't find how is this package build, my guess is that is as in the example of the OrganismDBI package. See that vignette section to build your own data package with the genome you want.

The locations are different because you ask for the transcript start and end, if you add the "TXNAME", you will (probably) see there are different transcripts, so different mappings are reported (untested).

Update

The example of the OrganismDBI package builds correctly the Homo.sapiens package, which appears to be linked with hg19. I provide here the relevant code to build the package for hg19 and hg38.

> source("https://bioconductor.org/biocLite.R")
> biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")
> biocLite("TxDb.Hsapiens.UCSC.hg38.knownGene")

> library(OrganismDbi)

# instructions for hg19, for hg38 change accordingly
> gd <- list(join1 = c(GO.db="GOID", org.Hs.eg.db="GO"),
             join2 = c(org.Hs.eg.db = "ENTREZID",
                       TxDb.Hsapiens.UCSC.hg19.knownGene = "GENEID"))
> destination <- tempfile()
> dir.create(destination)
> makeOrganismPackage(pkgname = "Homo.sapiens.hg19", graphData = gd,
                      organism = "Homo sapiens", version = "1.0.0",
                      maintainer = "Maintainer<maintainer@email>",
                      author = "Author Name", destDir = destination,
                      license = "Artistic-2.0")

Creating package in /var/folders/sf/tfrswpc96jl8b4520mn784440000gn/T//RtmpzR4YnY/file111b490af70f/Homo.sapiens.hg19

# installing package
> install.packages(destination, repos = NULL, type="source")

Then, you can load and use the packages with:

> library(Homo.sapiens.hg19)
> library(Homo.sapiens.hg38)

> select(Homo.sapiens.hg19, keys=c('KRAS'),
     cols = c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'),
     keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25358180 25403854
2   KRAS   chr12 25386768 25403863
> select(Homo.sapiens.hg38, keys = c('KRAS'), cols =
       c('SYMBOL', 'TXCHROM', 'TXSTART', 'TXEND'), keytype = 'SYMBOL')
'select()' returned 1:many mapping between keys and columns
  SYMBOL TXCHROM  TXSTART    TXEND
1   KRAS   chr12 25204789 25250931
2   KRAS   chr12 25209168 25250936
3   KRAS   chr12 25209431 25250803
4   KRAS   chr12 25233819 25250929


Is is possible to predict ncRNAs from sequence and homology alone? I'm working with a set of homologous genes (let's call it gene A) from several bacterial species.
I know (from previously published research) that in gene B (a close paralogue of gene A), there is a non-coding, antisense RNA. This antisense RNA regulates gene B by forming a duplex with the gene B mRNA.
I have noticed that, sometimes, an antisense RNA would be predicted for gene A in the same location in the IMG/MER database. This is not predicted for all of the available gene A sequences, only some of them. IMG/MER uses INFERNAL to predict RNAs, I think.
My question is - is there any way to assess how real this is without wet lab data? As it is, it could be that there's no actual antisense RNA in gene A, and this is a spurious prediction due to homology. Or it could be that the antisense RNA is present in all of the gene A sequences, and the RNA prediction software simply missed some of them. Is there a computational way of testing which of these options (if any) is more likely? Thank you for your time!

Is there any conserved secondary structure in this RNA? One could align all known instances of this RNA and see if there is any significant covariation according to R-scape. That could serve as a computational signal that the RNA is real (although wet lab experiments would still be needed).
Another question is whether this RNA is already in Rfam? If so, then the job of making an alignment of known instances is easier.


AND or OR gene - protein logical relation in KEGG How do I know if it is a AND or OR relationship between genes associated with the same enzyme in KEGG? E.g., for EC:1.6.1.2 (http://www.genome.jp/dbget-bin/www_bget?ec:1.6.1.2) does the bacterium Streptomyces coelicolor (SCO) need SCO7622 AND SCO7623 or SCO7622 OR SCO7623?

The enzyme EC:1.6.1.2 is an NAD(P) transhydrogenase and consists of two subunits. The two subunits are encoded by two different genes pntA and pntB.

So in terms of computer logic I would say AND, both subunits are necessary for a functional NAD(P) transhydrogenase enzyme (EC:1.6.1.2).


how to train a gene dataset with a nearest shrunken centroid classifier? I have a data file named "geneexp.csv".
the data contains information about gene expression of three different cell types (CD4 and CD8, CD19) I want to classify cells by performing the nearest shrunken centroid classification of training data in which the threshold is chosen by cross-validation. I split the data (70% train and 30% test).
data = read.csv("geneexp.csv")

splitData <- function(data, trainRate) {
  n <- dim(data)[1]
  idxs <- sample(1:n, floor(trainRate*n))
  train <- data[idxs,]
  test <- data[-idxs,]
  return (list(train = train, test = test))
}

split <- splitData(data, .7)
train <- split$train
test  <- split$test

then with the use of pamr package I tried to buid the following model and plot :
y <- train[[ncol(train)]]
x <- t(train[,-ncol(train)])
mydata <- list(
  x = x,
  y = as.factor(as.factor(y)),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

# Training and cross-validating threshold
model <- pamr.train(mydata)
cvmodel <- pamr.cv(model, mydata)
pamr.plotcv(cvmodel)

but I can't make it work. I get the following error:

Error in contrasts<-(*tmp*, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels

I have already transfer y to the factors. Can you help me? How can I fix it?

It's a bit different from usual R, but if you check the help:

data: The input data. A list with components: x- an expression genes
in the rows, samples in the columns)

So in your case, you need to transpose the matrix:
mydata <- list(
  x = t(x),
  y = as.factor(y),
  geneid = as.character(1:nrow(x)),
  genenames = rownames(x)
)

As I don't have your amazing data, I can only use iris below:
pamr.train(list(x=t(iris[,1:4]),y=iris[,5]))

pamr.train(data = list(x = t(iris[, 1:4]), y = iris[, 5]))
   threshold nonzero errors
1   0.000    4       6
2   0.841    4       7
3   1.682    4       10
4   2.523    4       11
5   3.364    4       13
6   4.205    4       18
7   5.046    3       22
8   5.887    3       23


What could be the reason for the samples not clustering? I'm performing RNA-seq analysis. I have used Hisat2 for aligning reads to the genome and stringtie for quantification and extracted read count information directly from the files generated by StringTie using "prepDE.py" mentioned in Stringtie manual. With the counts data I have used edgeR for Differential analysis.

Before that I have seen MDS plot showing in the following way:

I see that three samples are not clustered and thought that they are not sequenced well. So I removed the three samples and again made an MDS plot.

Can anyone tell me what could be the reason for this type of clustering? Is that because of sequencing? or using count data directly from stringtie or anything else?

Samples will only cluster by experimental group if the experimental effect is large enough that it's the primary source of variance between your samples. If that's not the case then you'll get results like you're seeing. Whether that's a problem or not ends up depending on how large of an effect you expect to see. If your groups are different stages in differentiation or another large effect, then you expect that samples will cluster into groups. If you're looking at something like some large unaffected tissue in a disease vs. normal state then seeing a cloud of points is quit reasonable.


tx2gene file for tximport issue I am trying to use tximport to get raw count matrices of a bulk RNA-seq sequencing experiment. I followed the tutorial:

https://bioconductor.org/packages/3.7/bioc/vignettes/tximport/inst/doc/tximport.html

Where they generate tx2gene in the following way:

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

So, I reasoned to do it the same way. But since my data is mouse one, I used a different gene annotation file:

library(TxDb.Mmusculus.UCSC.mm10.knownGene)
txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
k <- keys(txdb, keytype = "TXNAME")
tx2gene <- select(txdb, k, "GENEID", "TXNAME")

The transcript annotation file that I used when running salmon is gencode.vM13.transcripts.fa.sidx. Then when I try getting txi object:

txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

It raises the error:

Seems like I am creating tx2gene in the wrong way. Maybe I need to use a different library, not TxDb.Mmusculus.UCSC.mm10.knownGene. Any suggestions would be greatly appreciated.

  Update

I looked into the quant.sf files that I got from salmon and it looked the following way:

A colleague adviced me to strip away everything that follows the dot in the Name column, so it started to look:

However, still no luck. The identifiers in the k <- keys(txdb, keytype = "TXNAME") look like that:

"uc007afg.1" "uc007afh.1" "uc007afi.2" "uc011wht.1" "uc011whu.1"
"uc057aty.1"

I suppose I need to convert these identifiers first to ENS-like ones, but do not know how to.

The solution was to use a different library as I expected. It should be library(TxDb.Hsapiens.UCSC.hg19.ensGene). And also, as specified in the question, to remove everything after the dot in the Name column of the quant.sf files.


Nextflow: Runs without error, produce few but not all output files into output channel I was trying to develop a pipeline using nextflow and plink. My code runs without any error, but got only two output files instead of 22 (one for each chromosome). Here is part my code block, any help?
process step5 {
publishDir params.publish, mode: 'copy'
input:
val(hrc) from HRC_1000G_check_bim
val(hrc_site) from hrc_sites_file
tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
output:
tuple val("${prefix}"), val(chr_name), path("LOG-${prefix}${chr_name}.step4-HRC.txt") into step4_log
tuple val("${prefix}"), val(chr_name), path("Exclude-${prefix}${chr_name}.step4-HRC.txt") into step4_exclude
tuple val("${prefix}"), val(chr_name), path("Chromosome-${prefix}${chr_name}.step4-HRC.txt") into step4_chromosome
tuple val("${prefix}"), val(chr_name), path("Position-${prefix}${chr_name}.step4-HRC.txt") into step4_position
tuple val("${prefix}"), val(chr_name), path("Strand-Flip-${prefix}${chr_name}.step4-HRC.txt") into step4_strand
tuple val("${prefix}"), val(chr_name), path("Force-Allele1-${prefix}${chr_name}.step4-HRC.txt") into step4_allele
tuple val("${prefix}"), val(chr_name), path("ID-${prefix}${chr_name}.step4-HRC.txt") into step4_id
script:
"""
perl "${hrc}" -b "${prefix}${chr_name}.step4.bim" -f "${prefix}${chr_name}.step4.frq" -r "${hrc_site}" -h
"""
}

The output files are:
gsa-out/Chromosome-chr22.step4-HRC.txt     gsa-out/Force-Allele1-chr3.step4-HRC.txt  gsa-out/Position-chr22.step4-HRC.txt
gsa-out/Chromosome-chr3.step4-HRC.txt      gsa-out/ID-chr22.step4-HRC.txt            gsa-out/Position-chr3.step4-HRC.txt
gsa-out/Exclude-chr22.step4-HRC.txt        gsa-out/ID-chr3.step4-HRC.txt             gsa-out/Strand-Flip-chr22.step4-HRC.txt
gsa-out/Exclude-chr3.step4-HRC.txt         gsa-out/LOG-chr22.step4-HRC.txt           gsa-out/Strand-Flip-chr3.step4-HRC.txt
gsa-out/Force-Allele1-chr22.step4-HRC.txt  gsa-out/LOG-chr3.step4-HRC.txt

I am not sure why I am getting only 2, instead of 22 output files. With similar output channel previously generated 22 output files in the same script. For eaxmple:
process step1b {
publishDir params.publish, mode: 'copy'
input:
tuple val(prefix), val(chr_name), path(bfiles) from step1a_out
output:
tuple val("${prefix}"), val(chr_name), path("${prefix}${chr_name}.step1c.{bed,bim,fam}") into step1b_results
path("${prefix}${chr_name}.step1c.log") into step1b_logs
path("${prefix}${chr_name}.step1b.snplist.txt") into step1b_snplist
script:

"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${bfiles[1]
}" | grep ambig > "${prefix}${chr_name}.step1b.snplist.txt"
plink1.9 --bfile "${prefix}${chr_name}.step1b" --exclude "${prefix}${chr_name}.step1b.snplist.txt" --make-bed --out "${prefix}${chr_name}.step1c"
"""
}

I suspect your HRC_1000G_check_bim and hrc_sites_file channels are queue channels when they should be value channels. Please again see how multiple channels work. Also, third-party scripts, like your Perl script, would be better placed in a directory called bin/ at the root of your project repository. Nextflow adds this directory to your $PATH, so scripts in this directory will be automatically accessible in your pipeline without needing to localize them using a separate channel. This way also avoids having to specify an absolute path to invoke them. Just make sure your script is executable: chmod +x your_script.pl
hrc_sites = file( params.hrc_sites_file )

process step5 {

    input:
    tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
    path hrc_sites

    """
    your_script.pl \\
        -b "${prefix}${chr_name}.step4.bim" \\
        -f "${prefix}${chr_name}.step4.frq" \\
        -r "${hrc_sites}" \\
        -h
    """
}


Plotting gene models - alternatives to ggbio I am curious how people make genome track plots such as what can be done with the Bioconductor package ggbio.

For example this plot was made w/ ggbio:

I think ggbio is great for making fast plots, but it is kind of a pain to make publication quality figures because it is hard to manipulate the plots after they are constructed. Normally I would use ggplot2 to customize the plot extensively, but it doesn't work very well after plotting with ggbio.

How can I use ggbio for these types of plots? Is there a better tool instead?

EDIT:
I would get the legends for each individual plot and combine them before adding back to the main plot using something like this:

leg1.grob = cowplot::get_legend(plot1)
leg2.grob = cowplot::get_legend(plot2)
sharedleg.grob = gridExtra::gtable_rbind(leg1,leg2)

mainplot.grob = ggplotGrob(main_plot)
mainplot.grob = gtable_add_cols(mainplot.grob,
                                widths = sum(sharedleg.grob$widths))
mainplot.grob = gtable_add_grobs(mainplot.grob,
                                 grobs = sharedleg.grob,
                                 t=t,b=b,l=l,r=r)
cowplot::ggdraw(mainplot.grob)

It seems like the ggbio object already contains the corresponding ggplot object.
You should be able to extract the ggplot from the ggbio_obj@ggplot slot and extract your gtables.
Small example:
library(ggbio)
library(ggplot2)
library(patchwork)

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- as(c("chr17:41239394-41319151:+"), "GRanges")
p_cov <- autoplot(fl.bam, which = wh)
#> reading in as Bamfile
#> Parsing raw coverage...
#> Read GAlignments from BamFile...
#> extracting information...

df <-  data.frame(x= sample(1:100,size = 100,replace = TRUE))
p_hist <- ggplot(df) +
    geom_histogram(aes(x))

p_hist / p_cov
#> Error: Can't add `e2` to a ggplot object.
#> Backtrace:
#>     █
#>  1. └─patchwork:::`/.ggplot`(p_hist, p_cov)
#>  2.   └─GGally:::`+.gg`(e1, e2)
#>  3.     └─e1 %+% e2
#>  4.       └─ggplot2:::add_ggplot(e1, e2, e2name)
#>  5.         ├─ggplot2::ggplot_add(object, p, objectname)
#>  6.         └─ggplot2:::ggplot_add.default(object, p, objectname)

p_hist / p_cov@ggplot
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Created on 2021-02-05 by the reprex package (v1.0.0)
Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Red Hat Enterprise Linux
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Berlin
#>  date     2021-02-05
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version  date       lib source
#>  acepack                1.4.1    2016-10-29 [1] CRAN (R 3.6.0)
#>  AnnotationDbi          1.48.0   2019-10-29 [1] Bioconductor
#>  AnnotationFilter       1.10.0   2019-10-29 [1] Bioconductor
#>  askpass                1.1      2019-01-13 [1] CRAN (R 3.6.0)
#>  assertthat             0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
#>  backports              1.1.8    2020-06-17 [1] CRAN (R 3.6.0)
#>  base64enc              0.1-3    2015-07-28 [1] CRAN (R 3.6.0)
#>  Biobase                2.46.0   2019-10-29 [1] Bioconductor
#>  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor
#>  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor
#>  BiocManager            1.30.10  2019-11-16 [1] CRAN (R 3.6.0)
#>  BiocParallel           1.20.1   2019-12-21 [1] Bioconductor
#>  biomaRt                2.42.1   2020-03-26 [1] Bioconductor
#>  Biostrings             2.54.0   2019-10-29 [1] Bioconductor
#>  biovizBase             1.34.1   2019-12-04 [1] Bioconductor
#>  bit                    1.1-15.2 2020-02-10 [1] CRAN (R 3.6.0)
#>  bit64                  0.9-7    2017-05-08 [1] CRAN (R 3.6.0)
#>  bitops                 1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
#>  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.0)
#>  BSgenome               1.54.0   2019-10-29 [1] Bioconductor
#>  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.0)
#>  cli                    2.0.2    2020-02-28 [1] CRAN (R 3.6.0)
#>  cluster                2.0.8    2019-04-05 [2] CRAN (R 3.6.0)
#>  colorspace             1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
#>  crayon                 1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
#>  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.0)
#>  data.table             1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
#>  DBI                    1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
#>  dbplyr                 1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
#>  DelayedArray           0.12.2   2020-01-06 [1] Bioconductor
#>  dichromat              2.0-0    2013-01-24 [1] CRAN (R 3.6.0)
#>  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.0)
#>  dplyr                  0.8.5    2020-03-07 [1] CRAN (R 3.6.0)
#>  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.0)
#>  ensembldb              2.10.2   2019-11-20 [1] Bioconductor
#>  evaluate               0.14     2019-05-28 [1] CRAN (R 3.6.0)
#>  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.0)
#>  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.0)
#>  foreign                0.8-71   2018-07-20 [2] CRAN (R 3.6.0)
#>  Formula                1.2-3    2018-05-03 [1] CRAN (R 3.6.0)
#>  fs                     1.4.0    2020-03-31 [1] CRAN (R 3.6.0)
#>  GenomeInfoDb           1.22.1   2020-03-27 [1] Bioconductor
#>  GenomeInfoDbData       1.2.2    2020-04-01 [1] Bioconductor
#>  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor
#>  GenomicFeatures        1.38.2   2020-02-15 [1] Bioconductor
#>  GenomicRanges          1.38.0   2019-10-29 [1] Bioconductor
#>  GGally                 1.5.0    2020-03-25 [1] CRAN (R 3.6.0)
#>  ggbio                * 1.34.0   2019-10-29 [1] Bioconductor
#>  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.0)
#>  glue                   1.4.1    2020-05-13 [1] CRAN (R 3.6.0)
#>  graph                  1.64.0   2019-10-29 [1] Bioconductor
#>  gridExtra              2.3      2017-09-09 [1] CRAN (R 3.6.0)
#>  gtable                 0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
#>  highr                  0.8      2019-03-20 [1] CRAN (R 3.6.0)
#>  Hmisc                  4.4-0    2020-03-23 [1] CRAN (R 3.6.0)
#>  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.0)
#>  htmlTable              1.13.3   2019-12-04 [1] CRAN (R 3.6.0)
#>  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.0)
#>  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.0)
#>  httr                   1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
#>  IRanges                2.20.2   2020-01-13 [1] Bioconductor
#>  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.0)
#>  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.0)
#>  labeling               0.3      2014-08-23 [1] CRAN (R 3.6.0)
#>  lattice                0.20-38  2018-11-04 [2] CRAN (R 3.6.0)
#>  latticeExtra           0.6-29   2019-12-19 [1] CRAN (R 3.6.0)
#>  lazyeval               0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
#>  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.0)
#>  magrittr               1.5      2014-11-22 [1] CRAN (R 3.6.0)
#>  Matrix                 1.2-17   2019-03-22 [2] CRAN (R 3.6.0)
#>  matrixStats            0.56.0   2020-03-13 [1] CRAN (R 3.6.0)
#>  memoise                1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
#>  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.0)
#>  munsell                0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
#>  nnet                   7.3-12   2016-02-02 [2] CRAN (R 3.6.0)
#>  openssl                1.4.1    2019-07-18 [1] CRAN (R 3.6.0)
#>  OrganismDbi            1.28.0   2019-10-29 [1] Bioconductor
#>  patchwork            * 1.0.0    2019-12-01 [1] CRAN (R 3.6.0)
#>  pillar                 1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
#>  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
#>  plyr                   1.8.6    2020-03-03 [1] CRAN (R 3.6.0)
#>  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.0)
#>  progress               1.2.2    2019-05-16 [1] CRAN (R 3.6.0)
#>  ProtGenerics           1.18.0   2019-10-29 [1] Bioconductor
#>  purrr                  0.3.4    2020-04-17 [1] CRAN (R 3.6.0)
#>  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
#>  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.0)
#>  RBGL                   1.62.1   2019-10-30 [1] Bioconductor
#>  RColorBrewer           1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp                   1.0.4.6  2020-04-09 [1] CRAN (R 3.6.0)
#>  RCurl                  1.98-1.1 2020-01-19 [1] CRAN (R 3.6.0)
#>  reprex                 1.0.0    2021-01-27 [1] CRAN (R 3.6.0)
#>  reshape                0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
#>  reshape2               1.4.4    2020-04-09 [1] CRAN (R 3.6.0)
#>  rlang                  0.4.6    2020-05-02 [1] CRAN (R 3.6.0)
#>  rmarkdown              2.1      2020-01-20 [1] CRAN (R 3.6.0)
#>  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.0)
#>  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor
#>  RSQLite                2.2.0    2020-01-07 [1] CRAN (R 3.6.0)
#>  rstudioapi             0.11     2020-02-07 [1] CRAN (R 3.6.0)
#>  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor
#>  S4Vectors              0.24.3   2020-01-18 [1] Bioconductor
#>  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.0)
#>  sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi                1.4.6    2020-02-17 [1] CRAN (R 3.6.0)
#>  stringr                1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
#>  styler                 1.3.2    2020-02-23 [1] CRAN (R 3.6.0)
#>  SummarizedExperiment   1.16.1   2019-12-19 [1] Bioconductor
#>  survival               3.1-11   2020-03-07 [1] CRAN (R 3.6.0)
#>  tibble                 3.0.0    2020-03-30 [1] CRAN (R 3.6.0)
#>  tidyselect             1.0.0    2020-01-27 [1] CRAN (R 3.6.0)
#>  VariantAnnotation      1.32.0   2019-10-29 [1] Bioconductor
#>  vctrs                  0.2.4    2020-03-10 [1] CRAN (R 3.6.0)
#>  withr                  2.4.1    2021-01-26 [1] CRAN (R 3.6.0)
#>  xfun                   0.15     2020-06-21 [1] CRAN (R 3.6.0)
#>  XML                    3.99-0.3 2020-01-20 [1] CRAN (R 3.6.0)
#>  xml2                   1.3.0    2020-04-01 [1] CRAN (R 3.6.0)
#>  XVector                0.26.0   2019-10-29 [1] Bioconductor
#>  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.0)
#>  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor


How do I carry out an ancestry/admixture test on a single VCF file? This is a question from /u/beneficii9 on reddit. The original post can be found here.

Through the Personal Genome Project, I have had my whole genome sequenced by Veritas, and have it in the form of a single VCF file for the whole genome and one BAS file for each chromosome. The reference genome associated with the VCF file is hg19. It has been helpful in health data; for example, I discovered I'm homozygous for the non-functional variant CYP-2D6 gene (rs3892097), which can render several common medications useless, and helps explain why some medicines didn't really work for me. My doctor has found this information very helpful.
Unfortunately, I can't find any way of looking at admixture or ancestry. I've tried setting everything up using a combination of VCFTools, Plink1.9, and ADMIXTURE, but I can't get it to work. I think for ADMIXTURE you have to have a bunch of genomes sorted by geographical origin to compare your genome against, but I'm not sure how to do that, and what's online isn't very clear to me. So scratch that one off.
I've tried converting the file to 23andme format (and at this /u/psychosomaticism has been very helpful). I did that (though it seems there were problems because of the way the VCF file was set up). But the websites that take the data want you to point them to your 23andme account, and that doesn't really work if you only have the file. 23andme doesn't provide for people who had their whole genomes sequenced. They want you to give them a saliva sample like everyone else.
So, what can I do?

A modified implementation of Vivek's answer.

peddy is a Python package that samples an input .vcf at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.

I downloaded the .vcf and .vcf.tbi for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom .ped file, NA12878.ped, with the contents:

NA12878 HG001   0       0       2       0

At the command line:

$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped

The output files all have the prefix myvcf., here's myvcf.pca_check.png


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

Could I parse this output into --outfmt 6?

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


How can I find out if my gene of interest is duplicated? To identify homologous genes, I performed tBLASTn analyses against my genome of interest resulting in sequences of a closely related species. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome.

How can I determine if these gene copies correspond to true duplication events?

I aligned the newly found sequences to the reference genes and I found good alignments in more than one region (20% - 50% identity), but I am unsure if this is enough to consider them as duplicates.

Is there a threshold to help me determine if they are indeed duplications? How can I distinguish between a true duplication event and a pseudogene or a paralogous gene? The target genome belongs to a polyploid organism, although the genome assembly type is listed as haploid.

Any help will be greatly appreciated.

Thank you all very much for your responses, I am now beginning to understand how to find out if my genes are duplicated.

I performed some of the analyses that you suggested and I would very much appreciate your input on the results.

First, I followed @terdon's advice and I took the best hit for one of the query (source) genes in the target genome and I used this sequence as a query against the query (source) genome, i.e., I performed a reciprocal BLAST analysis (in my case it was a reciprocal BLASTp analysis). The first hit was my original query gene, which had an E-value of 0.0. For other sequences the reciprocal BLASTp hits were more variable, ranging from 3e-42 to 2e-154, but the first hit was always the query sequence. I also looked at the protein domains and they are conserved, i.e., both the query and the target sequences have the same domains. As for the exonic structure, it is fairly similar, for e.g., the query has 12 exons while the target sequence has 10. Exons are similar in length. So, in summary, when the potential duplicates from my target genome are blasted against the original source genome, the best hit is the original query sequence.

As for @zhanf's suggestion, I did a basic phylogenetic analysis using Clustal Omega (neighbour-joining tree without distance corrections) using four sequences: the original query sequence and the three gene copies (seq_5-1, seq_5-2 and seq_5-3) that I found in the target genome. This is the result:

I also compared each of the sequences against the original query:

seq_5-1 vs query: 34% identity (query length 780 aa, seq_5-1 length 400 aa)
seq_5-2 vs query: 42% identity (query length 780 aa, seq_5-2 length 580 aa)
seq_5-3 vs query: 40% identity (query length 780 aa, seq_5-3 length 490 aa)

So in summary, these sequences would appear to be paralogues, right?

Again, thank you very much for your help, I greatly appreciate it!

Leah

You might want to have a look at What is the difference between orthologs, paralogs and homologs?. Paralogs are duplicates by definition, and whether a gene is a pseudogene or still active has no bearing on whether it comes from a duplication.

In any case, to figure out what's going in in your target genome, you should use the target's genes. So, first get the best hit for each of your query genes in the target genome. Then, to investigate its evolutionary history, use that sequence (the one from the target genome) as  query against the target genome. That will let you assess how similar the target's hits are to one another and that's how you can attempt to figure out their evolutionary history.

Note that 20% identity is very low for closely related species. While there is no golden threshold, I would expect recent duplication events (and in closely related species, a duplication event that is in one genome and not the other will be recent) to produce much more similar sequences. In any case, you can't do this based on sequence identity alone.

Are the protein domains of the gene products conserved? Is the exonic structure of the gene conserved? You should also try taking the putative duplicate of geneA from genome2 and blast it against genome1 (your original source genome). Is the best hit the geneA from genome1, the query you used to find the putative duplicate in genome2? Or is there a better hit for it against another gene?

These are all questions to investigate and take into account when reconstructing the evolutionary history. You can't just look at identity and use some threshold to decide this sort of thing.

If you blast the potential dupicate from your target genome against your original source genome, is the best hit the duplicate?


Albacore basecalling running but outputs 0 reads I am trying to basecall data produced by the MinION using the SQK-LSK109 kit and FLO-MIN106 flowcell via the command-line. My version of albacore is the latest (v2.3.4). I tried running using the following command-line (the name of the sequencing run and fast5s are different and are default from machine output):

python3.6 /home/dkennetz/.local/bin/read_fast5_basecaller.py -i ./seqrun/fast5/ -t 1 -s ./good_reads/ -f FLO-MIN106 -k SQK-LSK109 -r -o fastq

-t 1 is just for this test. It is running and finishing almost instantly with 0 errors but outputting no reads. I have previously run this same command-line using the SQK-LSK108 kit and had no problems.

I then went on to validate the fast5s using the h5_validate tool:

python3.6 /home/dkennetz/.local/bin/h5_validate /home/dkennetz/.local/lib/python3.6/site-packages/h5_validator/schemas/single_read_fast5.yaml seqrun/fast5/0/seqrun.fast5 -v

The output from this is:

    Matcher <'UniqueGlobalKey', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

Error at /:

Matcher <'PreviousReadInfo', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

I read on the community page that this may be because the header is missing some relevant meta-data information? Any help on how to fix this would be greatly appreciated!

This is cross-posted on the community page for ONT here:

https://community.nanoporetech.com/posts/albacore-basecalling-runni

I believe you have multiread fast5 files, as generated by recent versions of MinKNOW. A conversion script is available to convert these to the older, one-read one-file format. AFAIK Albacore is being deprecated and no version will be made supporting the multiread fast5 files, as Albacore is being replaced by Guppy.

The command line to convert from multi read fast5s to single read fast5s is as follows:

$multi_to_single_fast5 -i /path/to/fast5/ -s /path/to/save/output/


Clustering with Seurat My data consists of two groups control and stim. I'm concerned that my clustering in the stim group is dominated by the unregulated genes. I can remove these genes and recluster but is there a way to then project the removed genes onto the new umap space so see if they are indeed all clustered together. They are no longer in the seurat object so I can't think of how to do this. I hope this makes sense.

Yes, this can be done, particularly if your concern is with the clustering / UMAP rather than the transformation process; I've done this recently with some Rhapsody single cell data that we generated. The main trick is adding 'return.only.var.genes = FALSE', as well as calculating a variable set that is a bit larger than the desired set.
Note that UMAP is a technique to place cells, rather than genes. If you have the information on where the cells should go, then plotting gene expression on top of that via FeaturePlot should work without any reprojection needed.
Here's what I did (with experiment-specific variables excluded):
sct.data <- SCTransform(sct.data, vars.to.regress="percent.mt",
                        variable.features.n = 3000 + length(genesToExclude),
                        return.only.var.genes = FALSE);

sct.data <- RunPCA(sct.data, verbose = TRUE,
                   features=setdiff(VariableFeatures(sct.data), genesToExclude));
sct.data <- RunUMAP(sct.data);
sct.data <- FindNeighbors(sct.data);
sct.data <- FindClusters(sct.data);

All transformed genes should be available for querying in the results.


What is the best way to address the question of doublets and multiplets in a single cell RNA seq data set? I have attached a histogram plot of the number of genes per cell in a single cell RNA seq data set of lung endothelial cells. I do not find a bimodal or multimodal distribution of the number of genes per cell within any particular condition. Doublets/multiplets should have an aberrantly high number of gene expression and if many, can lead to a bi or multimodal histogram plot of the number of genes detected per cell.  Can I conclude that I have an insignificant presence of doublets or multiplets in my data set from this plot and hence continue with downstream analyses?
I am reluctant to use automatic doublet-finding algorithms because they have a high possibility to remove cells in a transitional differentiating phase. I also think that the default threshold of 2500 genes per cell given by Seurat might sometimes lead to the loss of good cells depending on the cell type.

Apart from trajectory analyses, are there any other analyses on single cell RNA seq data that are sensitve to doublets/multiplets?
Thanks in advance for your help.

Indeed there's no satisfactory doublet/multiplet identification methods as of 2020.
The most reliable way to remove doublets is still manual removal of doublet subclusters guided by marker genes.
Not only trajectory analysis can be impacted by doublets, but clustering and embedding can also be affected for obvious reasons.


What are the labels in my rstats distance matrix? I have a distance matrix generated by hierfstat, thusly, with a link to the fasta file here:

library(adegenet)
library(hierfstat)
snps <- fasta2DNAbin('test.fa', chunkSize = 50)
gi <- DNAbin2genind(snps)
# manually define the populations for passing to pairwise.fst
p <- c('Botswana', 'Botswana', 'Botswana', 'Botswana', 'France', 'France', 'Vietnam', 'Vietnam', 'Uganda', 'Uganda', 'Uganda', 'Uganda', 'Vietnam', 'Vietnam', 'Laos', 'Laos', 'Laos', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam')
f <- pairwise.fst(gi, p, res.type = c('dist', 'matrix'))
as.matrix(f)

which produces

    1   2   3   4   5
1   0   0.2189008   0.225567409 0.1821518   0.259409722
2   0.2189008   0   0.130736953 0.1648034   0.191050772
3   0.2255674   0.130737    0   0.1669077   0.006396789
4   0.1821518   0.1648034   0.166907683 0   0.203931457
5   0.2594097   0.1910508   0.006396789 0.2039315   0

Is it safe to assume that 1 = Botswana, 2 = France, 3 = Vietnam, 4 = Uganda and 5 = Laos? i.e. that the distance matrix rows/columns follow the order that each population first appeared in p?

Is there a way to determine for certain, which row/column in the distance matrix corresponds to which population?

Disclaimer: I don't have adagenet and hierfstat installed, I just looked at the source code. It seems like pairwise.fst (the deprecated code from the adagenet package, which you seem to be using) already returns a matrix. Set res.type="matrix" and try calling rownames(f) and colnames(f) afterwards.


Can I use a regular liner regression model when I'm working with DNA methylation data? I'm new on working with genetics data and I'm just wondering if I can use the lm function when I'm building my models or do I have to use lmFit function in limma package to build the models and what is the difference between the two packages in terms of performing linear regression models

You can use either, but lmFit has the benefit of returning an object that can be used with eBayes() so you can pool information across genes/probes/whatever. lm() is a base R function applicable basically everywhere. lmFit() if from the limma package, so originally intended for microarray data, though these days pretty much everything omics is analyzed with it since sample sizes tend to be small so everyone want's to pool information for more statistical power.

BTW, if you have methylation data then you'll want to fit data in logit space.


Error given while trying to index a BAM file with Samtools Index - NO COOR? I am currently working on my own Metagenomic pipeline, utilizing Bowtie 2 to map. Bowtie 2 outputs a SAM file, which I convert to a .BAM and sort it using Samtools. When I try to utilize Samtools to index my .BAM file it gives me this error:

[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

Here is my code:

source /opt/asn/etc/asn-bash-profiles-special/modules.sh
module load samtools/1.3.1

cd ~/gz_files/soapdenovo
r=20
###SAM to BAM and sort

samtools faidx soapdenovo.fa
samtools view -b -T soapdenovo.fa -o sample${r}bowtiemapping.bam sample${r}bowtiemapping.sam
samtools sort -n -T ~/gz_files/soapdenovo/temp -o sample${r}bowtiemapping_sorte$
samtools index sample${r}bowtiemapping_sorted.bam

################LETS MAKE THE ABUNDANCE TABLES###########################$
samtools view -c -f 4  sample${r}bowtiemapping_sorted.bam > sample${r}alignment$

samtools idxstats sample${r}bowtiemapping_sorted.bam >  sample${r}idxstats.txt

#git clone https://github.com/metajinomics/mapping_tools.git
#python mapping_tools/get_count_table.py *.idxstats.txt > contig_counts.tsv
##less contig_counts.tsv

Any help would be appreciated!

Notes on pipeline

samtools import was replaced with samtools view
faidxcorrected
Instructions used are here on github
Suggestion for future trouble shooting is given here  by @winni2k

By now this will be of little interest to the original poster, but for future reference:
$ samtools sort -n -T ~/[etc] -o sample${r}bowtiemapping_sorted.bam [etc]
$ samtools index sample${r}bowtiemapping_sorted.bam
[E::hts_idx_push] NO_COOR reads not in a single block at the end 1652 -1.

This error message indicates (admittedly somewhat obscurely) that the reads are not sorted by coordinate — in particular, there are mapped reads appearing after the unplaced unmapped reads (those with RNAME *) that should be at the very end of the file.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.


Get gene annotation type from gene name I am having RNA-seq data set for C.Elegans. I do not know which annotation scientists used for naming the genes of it, but I need to search somehow for the ^ERCC and ^mt- genes like specified in the tutorial:

https://www.bioconductor.org/help/workflows/simpleSingleCell/

If I just implement grepl search for the character array of gene names that I have, I am always getting all values FALSE as if there are no spike-ins and mitochindrial genes.

Any suggestion would be greatly appreciated.

Sample gene name from the dataset: WBGene00023193

The tutorial you linked to does not mention ERCC or mt genes.

Before running grepl or any other search, you should first check your inputs. It looks like you have WormBase gene IDs if all your genes look like WBGene00023193. Therefore, your question is actually how to convert WormBase gene IDs to gene symbols.

Here is one potential solution: https://github.com/WormBase/website/issues/4297

After you convert, you have to make sure that you have ERCC spike-ins in your gene table (probably not) and you mitochondrial genes start with "mt-".


Where can I find gene expression data on one of the cell lines in NCI-60? I'm trying to retrieve the most expressed genes from a cell line from the NCI-60, lets say cell line NCI-H23, but I'm not sure where to find the gene expression data. My main goal is to create a network using that data and iRefIndex to be loaded into Cytoscape. Does anyone know where I can find the correct data or guide me in the right direction?

CellMiner portal has plenty of data about NCI-60 cell lines, inclusive microarray gene-expression.

Also, NCI-60 is a subset of the ~1000 cell lines covered by the GDSC project for which you have multiple types of measurements (e.g. copy-number, mutation, gene-expression).


How to identify unknow bacterium species from whole genome genetic sequence I am a Biochemist that is unfamiliar with bioinformatic tools and new to academia as a whole.
I am currently using ILLUMINA PE data, which I trimmed (Trimmomatic), corrected (Rcorrector) and assembled (SPAdes). I am now interested in using the genetic sequences from my contigs to identify the original organism (I know it is a bacterium).
I have tried using BLASTn (with default parameters) on one of my contigs. I let the query run for 30 minutes before stopping, nothing came of it.
What would be the best way to go about doing this?

You can use KmerFinder and Kraken2 for taxonomy assignment


How to apply function on bedfile based on a reference bedfile I have two bed-like files:

File1

chr1     1     10    .    1    Sample1
chr1    20    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    1    Sample3

File2 (reference)

chr1     1     5     .    0
chr1     5    20     .    1
chr1    50   120     .   -1

I'd like to find a way to apply a function (in my case a subtraction) between the first file and the second, using the second as a reference, and producing a results file containing also the indication of the sample name.
For example, in my case the results file would be:

chr1     1      5    .    1    Sample1
chr1     5     10    .    0    Sample1
chr1    20     50    .    1    Sample1
chr1    50    120    .    2    Sample1
chr1   120    180    .    1    Sample1
chr1     1      5    .    1    Sample2
chr1    70    100    .    2    Sample3

I was thinking about implementing it in pandas, but it would be long to do and probably not so efficient. I looked for bedmap but it seems to me that this operation can be carried on (with the sum operation) only within a single file during aggregation, not between two different bed files. So maybe it is not the right tool to use in this case.

Also, it would be important in the final file to retain the sample name.

You could munge your first file into a BED5 file:

$ awk -vOFS="\t" '{ print $1, $2, $3, $6, $5 }' /tmp/f1.bed | sort-bed - > /tmp/f1.bed5

Sort the second file, if necessary:

$ sort-bed /tmp/f2.unknownSort.bed5 > /tmp/f2.bed5

Take the union of the two sorted files:

$ bedops -u /tmp/f1.bed5 /tmp/f2.bed5 > /tmp/union.bed

Here's what the union looks like:

chr1    1       5       .       0
chr1    1       5       Sample2 1
chr1    1       10      Sample1 1
chr1    5       20      .       1
chr1    20      180     Sample1 1
chr1    50      120     .       -1
chr1    70      100     Sample3 1

We have one name scheme for reference element IDs (.) and a second name scheme for sample IDs (Sample*).

We will exploit this difference in name schemes when we apply operations in a certain order. For this example, we will subtract a reference element score from a sample element score.

We can partition the union, mapping the disjoint elements back to the unioned set:

$ bedops --partition /tmp/union.bed | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed
chr1    1   5   .;Sample2;Sample1   0.000000;1.000000;1.000000
chr1    5   10  Sample1;.   1.000000;1.000000
chr1    10  20  .   1.000000
chr1    20  50  Sample1 1.000000
chr1    50  70  Sample1;.   1.000000;-1.000000
chr1    70  100 Sample1;.;Sample3   1.000000;-1.000000;1.000000
chr1    100 120 Sample1;.   1.000000;-1.000000
chr1    120 180 Sample1 1.000000

Each of the fourth and fifth columns are semi-colon-delimited strings. Each element is ordered by the order in which it is discovered in a bedmap operation.

Both columns have the same ordering or same one-to-one correspondence of ID to signal score. This is very useful, because we can turn these into a pair of ordered arrays and apply some logic, depending on the answer we want back.

For example, here is some awk logic to parse those two strings and subtract a reference element score from a sample element score:

$ bedops --partition /tmp/union.bed \
    | bedmap --echo --echo-map-id --echo-map-score --delim '\t' - /tmp/union.bed \
    | awk -vFS="\t" -vOFS="\t" \
        '{ \
             i=split($4,a,";"); \
             j=split($5,b,";"); \
             if ((i==1) && (a[i]==".")) { pass; } \
             else { \
                 for (k=1; k<=i; k++) { \
                     if (a[k]==".") { refIdx=k; } \
                 } \
                 for (k=1; k<=i; k++) { \
                     if (k!=refIdx) { print $0, ".", b[k]-b[refIdx], a[k]; } \
                 } \
             } \
         }' \
    | cut -f1-3,6-8 \
    | sort-bed - \
    > answer.bed

The file answer.bed will look like this:

chr1    1   5   .   1   Sample1
chr1    1   5   .   1   Sample2
chr1    5   10  .   0   Sample1
chr1    20  50  .   1   Sample1
chr1    50  70  .   2   Sample1
chr1    70  100 .   2   Sample1
chr1    70  100 .   2   Sample3
chr1    100 120 .   2   Sample1
chr1    120 180 .   1   Sample1

This partitioning of intervals is different from the one in your question. Perhaps this is due to how you partitioned your set. In any case, perhaps this would give you an idea of how to apply custom operations over a disjoint partitioning.

The main idea is that this part of the awk script finds the score associated with the reference file ("File2"):

for (k=1; k<=i; k++) {
    if (a[k]==".") { refIdx=k; }
}

Any index other than refIdx corresponds to a Sample* element. That means we can do operations between Sample* scores and the specific reference score.

We also leave out mappings where a reference score does not associate with any Sample* element:

if ((i==1) && (a[i]==".")) { pass; }

Thus, all the operations are always between a pairing of a sample and the mapped reference element.

The cut and sort-bed statements at the end just clean up the result so that you get something that looks like what you are expecting at the end, which is also sorted correctly and ready to use for downstream set operations.


What is Allelic Imbalance Can anyone help me explain what allelic imbalance is, hopefully shortly?
Surprisingly, we cannot find any introduction online.

It might be used differently in different contexts, but generally speaking, in my world - allelic imbalance is when there's a difference in the level of gene expression from different alleles, usually through genetic (e.g. a variant in a promoter) or epigenetic mechanisms (e.g. one copy silenced, as in imprinted regions).

There are a few blog posts here with some more information, and plenty of papers (google "Allelic imbalance") if you want to dig deeper.


counting individuals with the present gene and absent gene for each population I have a folder containing 500 files from different individuals:

$ ls
HI.3577.005.Index_25.QPG4_KS2W-35_all_PA
HI.3577.005.Index_27.QPH4_SD1W-13_all_PA
HI.3577.005.Index_9.QPA2_KS2W-8_all_PA
HI.3577.007.Index_1.Q004B_SD1A-1_all_PA
HI.3577.007.Index_21.Q076C_SD2W-23_all_PA
HI.3577.007.Index_23.Q079_KS1A-15_all_PA

Each of these individuals belongs to a population (in the example above individuals belong to KS2W, SD1W, KS2W, SD1A, SD2W, and KS1A respectively). Each of my files looks as follows (this sample belongs to from "MK4" population):

$ head HI.3579.001.Index_9.Q089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 0
Ha16_00023451 0
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

or another sample from MK4 population

$ head HI.221.009.Index_9.P089_MK4-4R_all_PA
Ha8_00040788 1
Ha4_00024045 1
Ha4_00025366 1
Ha16_00022130 1
Ha16_00023451 1
Ha8_00040789 1
Ha4_00025367 1
Ha4_00024046 0
Ha16_00022131 1
Ha16_00023452 1

the first column is genes and the second columns refers to the presence of that gene (coded as 1) or absence of gene (coded as 0).

For all individuals belonging to each of the populations, I would like to count the total number individuals that have the gene and individuals that do not have the gene.

For example, for the 2 examples above, as they both belong to MK4 population my desired output is going to looke like this:

$ head MK4
Ha8_00040788   2  0
Ha4_00024045   2  0
Ha4_00025366   2  0
Ha16_00022130  1  1
Ha16_00023451  1  1
Ha8_00040789   2  0
Ha4_00025367   2  0
Ha4_00024046   0  2

The first column is the gene, the second column is the number of individuals in the MK4 population that have the gene and the third column is the number of individuals that do not have that gene.

I have made a file containing population IDs

$ head pop.txt
KS2W
SD1W
SD1A
KS1A

Can anyone suggest me how I could do this? thanks

I can show you how I would approach it in R. You'll have to adjust and customize it to your own data of course. I only had two files to play with, from your example. Especially for the merge step you should build a for loop around it so that it will merge every next file automatically, since you have 500 files.

Here an example for the KS2W files:

myFiles <- dir(pattern="KS2W")

myFiles
[1] "HI.3577.005.Index_25.QPG4_KS2W-35_all_PA.txt"
[2] "HI.3577.005.Index_9.QPA2_KS2W-8_all_PA.txt"

fileNames <- list()

for(i in 1:2){
    fileNames[[i]] <- read.table(myFiles[i], sep = " ")
}

merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

present <- rowSums(merged[,2:3])

zeros <- rowSums(merged[,2:3] == 0)

output <- cbind(present, zeros)

rownames(output) <- merged[,1]

output
              present zeros
Ha8_00040788        2     0
Ha4_00024045        2     0
Ha4_00025366        2     0
Ha16_00022130       1     1
Ha16_00023451       1     1
Ha8_00040789        2     0
Ha4_00025367        2     0
Ha4_00024046        0     2
Ha16_00022131       2     0
Ha16_00023452       2     0

To add more files with merge, you could add something like this:

merged <- merge(fileNames[[1]], fileNames[[2]], by.x = "V1", by.y = "V1", all = F, sort = F)

for(i in 3:30){
  merged <- merge(merged, fileNames[[i]], by.x = "V1", by.y = "V1", all = F, sort = F)
}

You'll have to adjust all the numbers and coordinates to fit your data.


Are gene names same across species? I have a bunch of gene names of Apis mellifera (specifically 194). I used these gene names as an input on STRING database to create a network for Drosophila melanogaster. 54 of those genes were also present in D. melanogaster and I received a second network. Is this a credible method of creating homologous gene network. Or should I use the sequences of the 194 genes and BLAST it to the geneset of the D. melanogaster and do it that way?

Gene names are not consistent across species. Most of the time genes are automatically annotated and the algorithms try to give names that are similar to their matches in other organisms, but not always. It's also more likely that genes will be automatically named/annotated based on their sequence, and manually named/annotated based on their function.

The further away you get from Homo sapiens & Mus musculus (and maybe Caenorhabditis elegans, Arabidopsis thaliana & Danio rerio), the less likely it is to find consistently named genes.