How do I convert an Excel file (.xlsx) to an .smi file format? How do I convert an Excel file(.xlsx) to .smi file format?
The original file I have is an excel file containing over 27,000 ligands. I am going to using the Padel Descriptor software to calculate the descriptors for the ligands however, Padel Descriptor just supports .smi files. Screenshot of the excel file added below for more clarity.

So a smi file is simply a text file formatted SMILES tab name newline.
Therefore what you need to do is

save for safety
make a new sheet
copy the SMILES column and the name column
save as tab separated file
in Finder/Nautilus/File Explorer rename the file extension to .smi


How to simulate nanopore reads? I have looked already here: Tools for simulating Oxford Nanopore reads
. This doesn't answer my question, because it lists a few Nanopore read simulators, but I have specific problems with each of them.

I am trying to simulate some Nanopore reads. I could use NanoSim, but that tool seems to be very difficult to use. I only have a reference genome, but you need training data as well? anyone have any experience with this? I can only find bacterial and yeast data (and I want to try this on Arabidopsis thaliana reference)

Another promising tool I found, was a paper for SNaReSim, reference here:
https://ieeexplore.ieee.org/document/8031171. However I cannot seem to find the tool itself, which is also mentioned in the other stackexchange post I mentioned above.

Does anyone have experience simulating PacBio reads with NanoSim, SNaReSim or maybe another tool? Because I am not sure about which one to use. Others like SiLiCO or ReadSim were updated 2 years ago and I believe that Nanopore technology has improved a lot in that time.

To clarify a bit on why I want to simulate data: I just want to test several tools about SV calling and simulating for a school project. All of the points mentioned in the comments I will use for the discussion, but I wanted to focus on 1 read simulator per technology and didn't know what to use, because Nanopore is difficult to simulate in that sense. That was already known to me.

It is important to train an error model as NanoSim does, as we do not fully understand the error processes involved in both the nanopore sequencing process and the basecalling process. Any sort of read simulator that does not use a species-specific model is just not going to produce realistic simulated reads.

As far as datasets go, I would recommend looking at this assembly and the associated data uploaded on EMBL:

  Michael, Todd P., et al. "High contiguity Arabidopsis thaliana
  genome assembly with a single nanopore flow cell." Nature
  communications 9.1 (2018): 541.

There are several data types there, as well as an assembly to train the models on.

Is it OK to train a model using reads from one individual and a reference assembly from another individual?

Generally, it is not good to map reads from one individual to the reference genome of another individual since this introduces something called reference bias. When we align reads to a reference genome generated from the same reads, or just from other reads from the same individual we remove this problem. All true reads will map to correctly assembled parts of the genome.

Why can we train a model for nanopore reads on a reference generated from ... nanopore reads? Won't we get a bad model for indels?

Yes, well, maybe. If the reference genome has been corrected to remove indels using a tool like pilon, then you can generally presume that the reference genome has very few to no indels present as a result of sequencing error from the Nanopore/PacBio reads. As a result, the model derived from aligning reads should have the correct error profile for indels.

Again, it is important that the long-read-assembly was corrected with Illumina shotgun reads from the same individual from which the long reads were derived. Otherwise the corrected reference would have a chimaera of indels and SNPs from multiple individuals.


How to get all PDB homologs from Uniprot (mapping + BLAST)? I'd like to create a dataset consisting of all sequences which are either present in the PDB, or whose homolog is present in the PDB. In other words, any sequence in the PDB or any sequence related to it. The similiarity margin is to be very wide, so anything above, for instance, 25% sequence identity is accepted (perhaps there are better criteria, but you get the point).

No further limitations, so any organism, with or without experimental data, not limited to the manually curated Swissprot set.

I have tried the SIFTS database, but that contains, as far as I can tell, only mappings for sequence identity >85%. I would like to go far below that.

For the correspondence between PDB and Uniprot entries you can use SIFTS -- a semi-automated mapping between PDB and UniProt maintained by PDBe.

The pipeline that creates the mappings uses BLAST and a few other criteria to decide which UniProt entry should be assigned to each PDB entry. The SIFTS website has all the data in CSV files.


BAM to BigWig without intermediary BedGraph I have a pipeline for generating a BigWig file from a BAM file:

BAM -> BedGraph -> BigWig

Which uses bedtools genomecov for the BAM -> BedGraph part and bedGraphToBigWig for the BedGraph -> BigWig part.

The use of bedGraphToBigWig to create the BigWig file requires a BedGraph file to reside on disk in uncompressed form as it performs seeks. This is problematic for large genomes and variable coverage BAM files when there are more step changes/lines in the BedGraph file. My BedGraph files are in the order of 50 Gbytes in size and all that IO for 10-20 BAM files seems unnecessary.

Are there any tools capable generate BigWig without having to use an uncompressed BedGraph file on disk? I'd like this conversion to hapen as quickly as possible.

I have tried the following tools, but they still create/use a BedGraph intermediary file:

deepTools

Some Benchmarks, Ignoring IO

Here are some timings I get for creating a BigWig file from a BAM file using 3 different pipelines. All files reside on a tmpfs i.e. in memory.

BEDTools and Kent Utils

This is the approach taken by most.

time $(bedtools genomecov -bg -ibam test.bam -split -scale 1.0 > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    1m20.015s
user    0m56.608s
sys     0m27.271s

SAMtools and Kent Utils

Replacing bedtools genomecov with samtools depth and a custom awk script (depth2bedgraph.awk) to output bedgraph format has a significant performance improvement:

time $(samtools depth -Q 1 --reference test.fasta test.bam \
  | mawk -f depth2bedgraph.awk \
  > test.bedgraph \
  && bedGraphToBigWig test.bedgraph test.fasta.chrom.sizes kent.bw \
  && rm test.bedgraph)

real    0m28.765s
user    0m44.999s
sys     0m1.166s

Although it is has less features, we used mawk here as it's faster than gawk (we don't need those extra features here).

Parallelising with xargs

If you want to have a BigWig file per chromosome, you can easily parallelise this across chromosome/reference sequences. We can use xargs to run 5 parallel BAM->BedGraph->BigWig pipelines, each using the tmpfs mounted /dev/shm for the intermediary BedGraph files.

cut -f1 test.fasta.chrom.sizes \
  | xargs -I{} -P 5  bash -c 'mkdir /dev/shm/${1} \
  && samtools depth -Q 1 --reference test.fasta -r "${1}" test.bam \
  | mawk -f scripts/depth2bedgraph.awk \
  > "/dev/shm/${1}/test.bam.bedgraph" \
  && mkdir "./${1}" \
  && bedGraphToBigWig \
    "/dev/shm/${1}/test.bam.bedgraph" \
    test.fasta.chrom.sizes \
    "./${1}/test.bam.bw" \
  && rm "/dev/shm/${1}/test.bam.bedgraph"' -- {}

deepTools

Let's see how deepTools performs.

time bamCoverage --numberOfProcessors max \
  --minMappingQuality 1 \
  --bam test.bam --binSize 1 --skipNonCoveredRegions \
  --outFileName deeptools.bw

real    0m40.077s
user    3m56.032s
sys     0m9.276s

This can be done in R very easily from an indexed .bam file.

Given single-end file for sample1.

library(GenomicAlignments)
library(rtracklayer)

## read in BAM file (use readGAlignmentPairs for paired-end files)
gr <- readGAlignments('sample1.bam')

## convert to coverages
gr.cov <- coverage(gr)

## export as bigWig
export.bw(gr.cov,'sample1.bigwig')

Be aware that this method doesn't include normalization steps (such as normalizing to total coverage). Most of these additional steps can be added if necessary.


How does FDRtool work? I have a question about using FDRtool. In the below code (on RNA seq data whose p values were acquired using Deseq2), the FDRtool was first used and thereafter p.adjust using the benjamini hochberg method. I didn't write this, its taken from someone else's code found online, in addition it is also the way the it is written in the pipeline at my university. Therefore, I assumed that this is how fdrtool is used and in my eyes it looks like they both adjust the p value with the same method. So that's why I asked, as I want to understand why you would need both these lines, which I've seen many people do:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = T)
ddsRes[,"padj"]  <- p.adjust(FDR.ddsRes$pval, method = "BH")

Dont they both correct for false discovery rate? Why do you need to run p.adjust after using fdrtool?
Thanks!!

Something like this book chapter might help on understanding FDR. The packages or functions you call differ in their estimation of pi0 or the proportion of null hypothesis.
Basically, the p.adjust(..method = "BH") is a more conservative method with Benjamini-Hochberg, you are assuming the proportion of null features to be 1.
When you use fdrtool, you are estimating the proportion of null hypothesis and from there, the local FDR. This in theory gives you more power.
However, I have to note that in your code, you are not using fdrtool to estimate the fdr as explained above, because:
FDR.ddsRes <- fdrtool(ddsRes$stat, statistic= "normal", plot = FALSE)

With this, the raw p-values are stored in FDR.ddsRes$pval and the q value, which is the corrected p-value, is stored in FDR.ddsRes$qval. In your code you are basically taking the raw p-values, and applying Benjamin-Hochberg:
p.adjust(FDR.ddsRes$pval, method = "BH")

You should apply a two tailed test. From the wald stats, it should be:
derived_p = 2*pnorm(-abs(ddsRes$stat))

In your code, you are running it wrongly by passing the test statistic into fdrtools, we can compare the raw p-values:
head(data.frame(stat = ddsRes$stat,derived = derived_p,
DESeq2=ddsRes$pvalue,fdrtools=FDR.ddsRes$pval))

        stat    derived     DESeq2  fdrtools
1  1.2405796 0.21476110 0.21476110 0.1216358
2  1.2314421 0.21815755 0.21815755 0.1244054
3 -0.3292516 0.74196552 0.74196552 0.6811997
4 -0.4964574 0.61957174 0.61957174 0.5356140
5 -0.4101163 0.68172066 0.68172066 0.6088427
6 -1.8007943 0.07173531 0.07173531 0.0246428

To summarize, use the adjusted p-value returned by DESeq2. The code you have above is wrong unless there's a good reason to do a one-tail test.


What are 2D reads in the Oxford MinIon? Reading the Mash: fast genome and metagenome distance estimation paper, I see the reference to "2d reads" notion under the minION context (a sequencing technique)
e.g

    In both cases Mash was able to
correctly differentiate these closely related species
(ANI ≈ 95 %) using 43,806 and 91,379 sequences collected
from single MinION R7.3 runs of B. anthracis
Ames and B. cereus ATCC 10987, respectively (combined
1D and 2D reads) <----------  HERE.

But couldn't find any reference in Wikipedia or online papers to what it isCan someone shed some light?

Early MinION sequencing runs had forward and reverse DNA templates joined together by a hairpin adapter, so that the sequencer would read both strands from the same template. The consensus sequence that was generated by combining the base calls of opposite strands is referred to as a 2D read.

Due to a dispute with Pacific Biosciences, and for other reasons related to chemistry and basecalling, the equivalent mode in current sequencing runs is called 1D². There is no physical ligation of forward and reverse templates, but the kit and flow cell chemistry are modified to encourage the reverse strand to follow the forward strand very soon afterwards through the pore, with the joining and consensus generation carried out in software after the sequencing is done. In the output fast5 files, these 1D² reads can also be referred to as 2D.


Variants in Exome that affect miRNAs related to Schizophrenia I have 3 exomes results (VCF files), I would like to study variants (SNP, INDEL) that can affect the miRNAs, which are related to a complex disease (Schizophrenia). What is the best way to assess the effects of variants on miRNAs and how can I link between the damaged miRNAs and Schizophrenia. Can I find a help or advice please?

I would run them through the Ensembl VEP, which will tell you which genes they hit and the effect they have on them. Then I would filter by consequence "mature_miRNA_variant", which will find you all the variants in your file that fall in the mature miRNAs. Or you could filter instead by the gene biotype "miRNA" which will tell you if the variants fall anywhere in the miRNAs (which will include the full stem loop structure). I would then export the list of miRNAs and try to find out about their targets using miRbase.


Why do I have >10,000 cells in the 10X matrix produced by cellranger? In my scRNA seq experiment, single-cell libraries were generated using the GemCode Single-Cell Instrument and Single Cell 3′ Library & Gel Bead Kit v2 and Chip Kit (10x Genomics) according to the manufacturer’s protocol. About 8700 cells were added to each channel with a targeted cell recovery estimate of 5000 cells.

However, after processing my raw data with CellRanger software, I get 737,280 cells (barcodes) and about 31,000 genes. I cannot understand why I get so many cells and why this barcode whitelist is not up to 10,000 cells. After filtering this data retaining genes expressed in >= 5 cells as well as retaining cells with at least 200 detected genes, I get around 2200 cells.

Thanks

The raw data from cell ranger contains all of the barcodes detected in the experiment.  These raw barcodes are filtered by cell ranger to identify the barcodes that likely represent cells rather than empty droplets / dead cells. You can find the filtered data in the filtered feature-barcode matrices.  You can also examine the Run summary HTML file to see a summary of the cell ranger barcode filtering.

The reason why there are so many raw barcodes is an important optimization by 10X to increase the cell capture rate compared to techniques such as drop-seq.  In a 10X experiment the number of gel beads is much higher than the number of cells which leads to a higher probability of a cell and a bead ending up in the same droplet.  This also leads to a large number of empty droplets with just a bead and some ambient RNA or a dead cell.  These all end up being sequenced and the empty drops tend to have a small number of unique molecules (the majority less than 10 UMIs / barcode).


Demultiplex nanopore reads with custom barcodes We have a problem trying to demultiplex MinION sequences with custom barcodes. Do you have any software recommendations we can try for demultiplexing or how to demultiplex these custom barcodes with Albacore? We have tried using albacore but it only recognizes barcodes provided by Nanopore Technologies.

My semi-automatic workflow for demultiplexing nanopore reads using LAST can be found here - pay attention to the comments, where I've identified bits in the steps which might trip up people trying to follow them:

https://dx.doi.org/10.17504/protocols.io.7vmhn46

While this protocol is designed around using the ONT barcodes, the barcode fasta file can be replaced by another file [using barcode ID in the header] without changing the protocol. Make sure that this file only includes the unique parts of the barcodes (with a few bases of context if necessary to exclude random matches), so that reads won't be matched to the non-unique part of the wrong barcode sequence.


how to interpret of "pvclust" dendrogram and finding height for cutting dendrogram? I study on RNA-seq expression dataset about one cancer in TCGA. I downloaded FPKM dataset and removed batch effect by ComBat() function. Now, I used pvclust for clustering 72 samples and the result of that can access in this link.my code for running pvclust is as below:

fit_LMS <-pvclust(datForpvclust,method.hclust = "average",
                          method.dist = "euclidean",nboot=10000, parallel=TRUE )
pvrect(fit_LMS, alpha=0.95)

As you see, by alpha=0.95 just 3 clusters are significant although in right cluster a few clusters have been included. really I can't analysis this result.
So, I have 2 questions:
1-  What is the interpretation of pvclust dendrogram? Does my dataset has meaning full clusters?

2-  I am interested in the height of tree cut in the dendrogram, which height is better for this dataset based on pvclust result? H=105 or H=110 or another height?

I appreciate it if anybody shares his/her comment with me.

You do not need to give links to your images, you can embed images in your questions using the "image" icon.
FPKM (and RPKM) might not be the best RNA-seq normalization method, see this and this.

pvclust providing p-values for the clusters computed by hclust. To do so, bootstrap resampling is used. From their documentation:

  Values on the edges of the clustering are p-values (%). Red values are
  AU p-values, and green values are BP values. Clusters with AU larger
  than 95% are highlighted by rectangles, which are strongly supported
  by data.

Basically you are getting a clustering pattern and a statistic associated with each of the clusters giving you an idea about the "statistical significance" of the clusters you have obtained (i.e. clusters most probably not being a result of strong batch effetcs, ...).

To be able to comment on the resulting clusters making sense or not, you will need to consult somebody with domain knowledge, it is not a number or two that would determine the quality of the output but also how good is the fit between your result and the (expected) biology. For example in your case with TCGA data:

Do the samples cluster by batch?
Do tumor and adjacent normal tissue samples aggregate into different clusters?
Do the samples in the same cluster share the same mutation/methylation pattern/disease subgroup?


What do the FASTQ file names mean here? I would like to run an analysis on Nanopore data. I'm trying to download sample data from https://github.com/nanopore-wgs-consortium/chm13. I would like to run minimap on the FASTQ files, resuling in a single BAM file.

Q: What do the file names mean? What is "S13"? What is "L002_l1"? What is "L002_R1"? Am I supposed to download them and combine them into a single FASTQ file, and then align them?

What is "S13"?

S# is most likely the sample number. If so, you do not want to combine these.

What is "L002_l1"?

L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).

What is "L002_R1"?

R1 is one end of paired-end sequencing, R2 is the other end. These filepairs must have the same number of lines. Use an aligner that utilizes both.


Is it wise to use RepeatMasker on prokaryotes? I'm looking for a way to identify low complexity regions and other repeats in the genome of Escherichia coli. I found that RepeatMasker may be used for example when drafting genomes of prokaryotes (E. coli example). But RepeatMasker works on a limited dataset of species, neither of them being prokaryotes. By default, when running RepeatMasker, if no species is specified, it will compare with homo sapiens data.

This seems rather inadequate, but the most relevent alternative, PRAP, requires a "dead" tool (VisCoSe, by Michael Spitzer).

Is it still wise to to use RepeatMasker on Escherichia coli?
If yes, which settings would maximise relevance ?

If I understood correctly your question, you want to mask those regions in a (FASTA?) genome. I think you could identify those regions using mummer and mask them using bedtools.
# align genome against itself
nucmer --maxmatch --nosimplify genome.fasta genome.fasta

# select repeats and convert the corrdinates to bed format
show-coords -r -T -H out.delta | awk '{if ($1 != $3 && $2 != $4) print $0}' | awk '{print $8"\t"$1"\t"$2}' > repeats.bed

# mask those bases with bedtools
bedtools maskfasta -fi genome.fasta -bed repeats.bed -fo masked.fasta

Have a look at nucmer and bedtools maskfasta options to fine-tune your analysis.


Perplexity is too large I am trying to run {Rtsne}. On pca dataset prepared by using dudi.pca from {ade4}.

For pca calculation I selected 40 components while specifying parameters for dudi.pca:

pca = dudi.pca(counts(sce_normalized), scannf=FALSE, nf=40)

{Rtsne} gives me an error: “perplexity is too large”, how could it be fixed?

The code I am trying to run is below:

name_of_file <- "dudi_pca.rds"
path_to_file <- paste0("/scratch/nv4e/scRNA-
seq/normalization/new_datasets/",
                   name_of_file, sep="")
dudi_pca <- readRDS(path_to_file)

# run Rtsne with default parameters
rtsne_out <- Rtsne(as.matrix(dudi_pca), pca=FALSE)

I was using dudi_pca incorrectly. The supplied parameter to as.matrix() should have been dudi_pca$li:

rtsne_out <- Rtsne(as.matrix(dudi_pca$li), pca=FALSE)


polishing assembled genome to QV50 value I am trying to assemble some genomes from pacbio reads using HGAP3 followed by quiver polishing. However, even after multiple rounds of quiver polishing, I cannot attain quality value 50 (QV50) for a genome.

QV50 seems to be recommended and I am able to get upto QV30 only. I have done polishing for almost 15 times already. I have tried polishing with different parameter changes like subread length cutoff and blasr maximum divergence.

Can anyone suggest me on how I can achieve an assembled genome of QV50 value?

Thank you

15 times is too many times. If you were going to get quiver to get you good QVs it would have happened after 2-3 rounds probably.

Some questions:

Why are you using quiver? Are your reads from RSII chemistry? Arrow is recommended for Sequel reads.
What is your dataset coverage? If you do not have very high coverage (I think >30X is usually a minimum for a eukaryotic genome) it may be that you are just not going to get exactly the values you want because assembling low coverage data is hard.
Have you tried other polishing tools? They might be able to get you better value. I am not an expert but some people like to use Racon or Apollo.
While I'm not an expert, looking into how the QVs are calculated can be helpful. I looked at the quiver docs and it looks like if the reads map ambiguously then QV is low. So if you for example have lots of haplotigs, then mapping might be ambiguous for a lot of reads, and therefore this will depress your QV values. This is therefore a problem that can come from the starting assembly and not just the polishing. You may want to go back and take another look at your starting assembly.


Run Nextflow with file dependencies inside a Docker container? I’m trying to run a Nextflow workflow in a custom Docker container. Without Docker, the workflow succeeds. But running it inside the container leads to an error because a dependent file cannot be found.

Here’s a minimal script to illustrate the issue (similar to the basic Nextflow script example):

params.reference = "/data/reference/GRCh38.fa"

reference = file(params.reference)
fastq_channel = Channel.fromFilePairs('/data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "$reference" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t "${reference}.bwt"
    """
}

Running this inside a container will cause the path of the query, /some/fastq/, to be mounted inside the container. However, the path to the database, /data/reference, will not be mounted. Consequently, I get the following execution error:

  stat: cannot stat 'GRCh38-chr21.fa.bwt': No such file or directory

Without the container, the command succeeds.

I understand why the error occurs: Nextflow doesn’t track the dependency on $reference inside the command, and thus doesn’t mount its path in the container.

My question is, how do I fix this? The documentation is silent on this regard, and none of the workflows I can find online address the problem. For example, the GATK4 IARC workflows would presumably have the same issue, and aren’t configured for use with Docker.

Just adding the reference as a process input is insufficient, since BWA requires various other index files. Am I supposed to add these all manually as inputs?

The following workaround seems to work but ideally I don’t want to have to change the script portion of my process so that it works with Docker.

params.ref = 'data/reference/GRCh38.fa'

ref = file(params.ref)
ref_index_dir = ref.parent
fastq_channel = Channel.fromFilePairs('data/fastq/*_R{1,2}.fastq.gz')

process map_reads {
    input:
    set sample, '?.fastq.gz' from fastq_channel
    file ref
    file ref_index_dir

    output:
    file("${sample}.sam") into sam_channel

    script:
    """
    echo 'bwa mem -Y "reference/$ref" 1.fastq.gz 2.fastq.gz >${sample}.sam'
    stat -t reference/${ref}.bwt && touch ${sample}.sam
    """
}

Here, reference inside the script section refers to the parent directory name of the reference, and needs to be adjusted accordingly. Having to refer to reference/$ref is annoying but at least it works — both inside a Docker container and outside of it.


Display a score or graph to visualize the degree of conservation of each residue from alignment data I would like to display a score or graph to visualize the degree of conservation of each residue from the amino acid alignment data.
If possible, I'd like to extract the parts of the game that scored above a certain level.
Is there a tool or source code (python) to do that?
Sorry for my bad English ...

I would suggest a logo plot, e.g.:

I would take @Devon Ryan's suggestion to create a multiple sequence alignment (MSA), and plug that into that logo plot tool (or others, there are many like it).


What do 0s in stageR mean? I'm running stageR as part of a pipeline based on DRIMSeq and wanted to know what it means when the transcript level padj is equal to 0, particularly when the gene is non-significant. I'm pretty sure that in some cases if the padj is 0 it means that the chance of it being different due to chance is 0. There are a lot more 0s than I would expect in sites where the gene padj is non-significant. One example, I have a gene padj that's >0.9 but both of the sites I'm looking at have a padj of 0.

I have a couple of suspicions about why this is: either because it defaults to this if the comparative difference is or it uses 0 as a substitute to NA. Am I thinking along the right lines? Are the 0s real or are they likely to be a quirk of the package?

In stageR, we automatically return stage-wise adjusted p-values of zero for transcripts coming from genes that were found to be significant and only have 1 transcript, when testing for differential transcript expression (DTE).
Similarly, when testing for differential transcript usage (DTU), we return p-values of zero for transcripts derived from a significant gene consisting of no more than two transcripts (based on the tx2gene object you provide).

Transcript-level p-values for non-significant genes should not occur. We provide NA's for genes with only one transcript for DTU, and by default you should get NA p-values for non-significant genes, since these genes do not proceed further than the screening stage.

Could you provide some more details on the analysis? Which multiple testing correction have you selected?

Would it be possible to provide a reproducible example for your analysis?


How can I download the SNP genotyping information from GTEX? I am new to the field of Bioinformatics. GTEx database is great; I know how to check whether a SNP is eQtl with the online tools of GTEx.

But my question is how to download the right data to draw a boxplot on my own. like this paper:

Please help me . I have been stuck on this for many days.

There are two types of data you can get from GTEx. The first is summary statistics, basically the slopes and p-values you see in the boxplot. You can also get expression levels for all genes measured in GTEx.

However, in order to plot boxplot by yourself, you also need genotype data. This is not publicly available because of privacy issues. In order to get it, you have to apply to it on dbGaP. I think you can find relevant information on the GTEx website.


Same transcript coordinates in gtf file, different transcript ID I have a gtf file from Ensembl, and I noticed that several "transcript" annotations have the exact same coordinates. See for instance the third and fourth transcripts ("Y74C9A.2b.1" and "Y74C9A.2b.4") for this gene:

$ grep "WBGene00022276" genes.gtf | awk '$3 == "transcript" {print}'
I   ensembl transcript  10413   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.2"; transcript_name "Y74C9A.2a.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS7921";
I   ensembl transcript  11495   16793   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.2"; transcript_name "Y74C9A.2b.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9220";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.1"; transcript_name "Y74C9A.2b.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.4"; transcript_name "Y74C9A.2b.4"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11505   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.3"; transcript_name "Y74C9A.2a.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS16820";
I   ensembl transcript  11618   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.1"; transcript_name "Y74C9A.2a.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS43845";
I   ensembl transcript  11623   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.3"; transcript_name "Y74C9A.2b.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS49800";

What does it mean? Do these two "transcripts" actually correspond to an exon, belonging to two different isoforms?

As Ian explained, these are different transcripts which happen to have the same start and end positions. You have no information on their exonic structure in that file. However, if you look them up at EnsEMBL, you will see:

Transcript Y74C9A.2b.1 :

1   Y74C9A.2b.1.e1  11,499  11,561  -   -   63
    Intron 1-2      11,562  11,617          56
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

This has 4 exons, two of which are non-coding (UTR). The second starts at

Transcript Y74C9A.2b.4

1   Y74C9A.2b.4.e1  11,499  11,557  -   -   59
    Intron 1-2      11,558  11,617          60
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

Granted, the two transcripts are very simlar, they both have 4 exons, the first 2 of which are non coding (UTR). However, the 1st exon of transcript
Y74C9A.2b.1 is 63 nt long and ends at 11,561 while the 1st exon of transcript Y74C9A.2b.4 is 59 nts and ends at 11,557.

So there you go: two (slightly) different transcripts.


Counting the number of paralogues for mouse genes gives me the wrong frequency in R I am trying to count the number of paralogues for the mouse homologues of the human protein-coding genes using BioMart. But for example in the 'PLIN4' gene its counting 35,000 paralogues instead of 4.

We think it is because some genes have one to many paralogues which causes repeats. When I run a single gene its gives me back the correct number of paralogues. Is there a way to either remove these repeats from the results or a way around this so that BioMart doesn't output these repeats.

I have also thought of maybe running one gene at a time, then counting it by setting up some sort of loop so that it does all of the genes from the list automatically.

The code I have written so far is:

# Load the biomaRt package:

library(biomaRt)
ensembl_hsapiens <- useMart("ensembl",
                          dataset = "hsapiens_gene_ensembl")
ensembl_mouse <- useMart("ensembl",
                       dataset = "mmusculus_gene_ensembl")

# Get all human protein coding genes:

hsapien_PC_genes <- getBM(attributes = c("ensembl_gene_id",
                                         "external_gene_name"),
                          filters = "biotype",
                          values = "protein_coding",
                          mart = ensembl_hsapiens)

ensembl_gene_ID <- hsapien_PC_genes$ensembl_gene_id

# Get mouse homologues

mouse_homologues <- getBM(attributes = c("ensembl_gene_id", "external_gene_name",
                                       "mmusculus_homolog_associated_gene_name"),
                        filters = "ensembl_gene_id",
                        values = c(ensembl_gene_ID),
                        mart = ensembl_hsapiens)

# Get mouse external gene name

mouse_homologues_external_gene_names <- mouse_homologues$mmusculus_homolog_associated_gene_name

mouse_paralogues <- getBM(attributes = c("hsapiens_homolog_associated_gene_name",
                                       "external_gene_name",
                                       "mmusculus_paralog_associated_gene_name"),
                        filters = "external_gene_name",
                        values = c(mouse_homologues_external_gene_names) , mart = ensembl_mouse)

# Remove genes with no paralogues
mouse_paralogs_data <- mouse_paralogues[!(is.na(mouse_paralogues$mmusculus_paralog_associated_gene_name)
                                          |
mouse_paralogues$mmusculus_paralog_associated_gene_name==""), ]

# Count paralogues per gene

library(plyr)
count_mouse_paralogues <- count(mouse_paralogs_data, "external_gene_name")
View(count_mouse_paralogues)

I'd use dplyr for this.

First, mouse paralogues contains no values for mmusculus_paralog_associated_gene_name that are NA:

library(dplyr)
mouse_paralogues %>%
  filter(is.na(mmusculus_paralog_associated_gene_name)) %>%
  nrow()

[1] 0

But there are values which are empty strings, so you can filter for those if you wish.

count_mouse_paralogues <- mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name != "") %>%
  count(external_gene_name)

Confirm that this gives the expected result (n = 4) for Plin4:

count_mouse_paralogues %>%
  filter(external_gene_name == "Plin4")

# A tibble: 1 x 2
  hsapiens_homolog_associated_gene_name     n
  <chr>                                 <int>
1 Plin4                                     4

And check mouse_paralogues:

mouse_paralogues %>%
  filter(hsapiens_homolog_associated_gene_name == "PLIN4")

  hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
1                                 PLIN4              Plin4                                  Plin3
2                                 PLIN4              Plin4                                  Plin5
3                                 PLIN4              Plin4                                  Plin2
4                                 PLIN4              Plin4                                  Plin1


Coverage calculation: long reads (RNA-seq) Say your aim is to calculate the coverage of an RNA-seq experiment generated with long-read sequencing (so, uneven read length).

Up to now, I relied on the Lander/Waterman equation:

$$C = L*N / G$$

where:

$C$ = final coverage
$G$ = haploid genome (or transcriptome) length
$L$ = read length
$N$ = number of reads

I have two major conceptual issues with this formula:

$L$ is uneven in this case. Should we use the median/mean read length instead? If so, side question: if the RNA-seq experiment was executed following the size-fractionated protocol (PacBio), is it advisable to merge all the fractions prior to coverage calculation or should it be calculated separately for each fraction?
Should $N$ consist of the mapped reads to the reference or simply the number of reads as found in the subreads whole dataset?

The read length is irrelevant when calculating the mean coverage statistic. It's simply the total number of bases sequenced divided by the target Xome length.

In the example provided in the question, $L*N$ is equivalently expressed as $\sum^{N} l_x$, or the total length of sequenced reads. For equal-length reads, this is easier to calculate as number of reads $*$ length of reads, but the full sum is necessary if read lengths are not equal.

For what it's worth, this value, $L*N$, is also equivalent to the average read length multiplied by the number of reads, i.e.

$$\frac{\sum^{N} l_x}{N} * N$$

... which makes it a bit more obvious why the total length of sequenced bases is very slightly easier to calculate.

In answer to the second part of the question, $N$ should ideally correspond to only the sequences that are from the reference. For an assembly, only the reads that were used to create the assembly should count when considering coverage. When there's no likelihood of contamination (and especially for RNASeq with variable coverage per transcript), a statistic using the total number of reads should be good enough for most purposes that I can think of. High-throughput sequencing is such a fuzzy process that a few tens of millions of bases here or there shouldn't make a difference.


mapping heteryzygous kmers on a genome I have a set (a couple of millions) of kmer pairs of known length (usually 21) that I know that are heterozygous in the middle nucleotide. For instance:

AAAAAAAAAANAAAAAAAAAA

Where N has always two states. I would like to map these nucleotides on a fragmented genome assembly that is presumably haploid. I bet that this problem is resolved, but I have hard times to make a practical solution, preferably within python environment. I considered following

Using an existing mapper

I tried to use python binding of minimap2, thinking that I just need to tell it that I search for exact matches only, but I did not manage to set it up. I set up smaller minimisers but it always returned no hits for kmers of length 21 (I tried to tweak parameters k, w and min_cnt, but to be honest I don't have a good understanding of how minimizers work). Is there a mapper that would be easily set up to do these searches?

Hash table

The second thing I tried was to build a hash table of all k - 1 / 2 mers of the genome to map the left and right sides of the kmers. Where sequences were mapped to lists of corresponding genomic positions, but this have blown my memory very soon (~50m of kmers were indexed using ~8GB of memory; i.e. this approach does not seems to be viable on a genomic scale).

I am also suspecting that there will be an elegant solution using a suffix array, but I have not cracked that one yet.

I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking bwa mem to report only perfect matches. So in your case, the idea is to make two copies of each 21-mer, replace the N with the appropriate allele, and then map the resulting sequences to your bwa-indexed assembly like so.

bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta

The -k 21 says to use a minimum seed length of 21, which forces exact matches. The -T 21 requires a minimum score of 21, which also enforces an exact match. The -a parameter reports all matches, since a "best" match doesn't make sense in this situation. The -c parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.

UPDATE If you don't want to introduce a dependency like bwapy, it's fairly straightforward to implement this in pure Python.

with tempfile.TemporaryFile() as kmerfile:
    for kmer in kmers:
        print('>{k}\n{k}'.format(k=kmer), file=kmerfile)
    kmerfile.flush()             # Often works without these two lines
    os.fsync(kmerfile.fileno())  # But best to be safe
    bwacmd = ['bwa', 'mem', '-k', ksize, '-t', ksize, '-a', '-c', '5000', refrfile, kmerfile.name]
    bwaproc = subprocess.Popen(bwacmd, stdout=subprocess.PIPE, universal_newlines=True)
    for line in bwaproc.stdout:
        # process SAM records line by line

Also, one can call bwa mem via a C API, but the mechanics of this are left as an exercise for the reader.


How can I use Arlequin via the command line? I've got a decent knowledge of programming (incl. bash scripting) but I fail to understand how Arlequin works. Could you please help me with a very simple reproducible example on how to use Arlequin via the command line?

As the .zip file comes with its pre-compiled version, I'll just note that I'm on Mac OS X but of course, a reproducible example with Linux (or at worst with Windows) would also be helpful.

For MAC OSX, the executable found in arlsumstat_macosx (arlsumstatmac_64bit) seem to work appropriately. The directory for arlsumstat_macosx does not contain any examples but you can use example files from arlecore_macosx/Example files_linux using this syntax.

The complication is that arlsumstatmac_64bit will look in the current directory for the file arl_run.ars. There is a file with the same name in arlecore_macosx but it is not compatible. So you have to be in the arlsumstat_macosx.

In short, cd to arlsumstat_macosx and do

./arlsumstatmac_64bit path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.ars path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.arp

, where the .arp file contains the data and the .ars file contains the settings (the description of what you want Arlequin to do).


Rejected publication A publication that we have been working on for over 6 months on COVID-19 has been rejected by 2 journals. The objective of our research was to identify effective compounds against the COVID-19 virus by using Molecular Docking studies performed with AutoDock Vina. Which Journal should we submit our manuscript to? Following are the Rejection emails I received from them. What do you think is the reason for the rejection and what should I do to resolve the issue to get my manuscript successfully published.
Rejection Notification from International Journal of Molecular Sciences
"We are writing to inform you that we will not be able to process your
submission further. Submissions sent for peer-review are selected on the
basis of discipline, novelty and general significance, in addition to the
usual criteria for publication in scholarly journals. Therefore, our decision
is not necessarily a reflection of the quality of your work."
Rejection Notification from International Journal of Biological Macromolecules
"Thank you for submitting your manuscript to International Journal of Biological Macromolecules.
I have reviewed your manuscript and I must regretfully inform you that your work presents information that is specific to a selective audience in its nature and interest and is more appropriate for a more specialized journal. The focus of your present work, a computational study, is unfortunately outside the scope of the journal, therefore, unsuitable for publication in IJBIOMAC.
For IJBIOMAC publication, the computational analysis should be used only as a supporting tool to non-computational analysis.
We appreciate you submitting your manuscript to International Journal of Biological Macromolecules and thank you for giving us the opportunity to consider your work."
Thank you for your time.

I cannot comment on your specific manuscript since I obviously cannot read it but in general the question is whether a study does add an interesting novelty.
The problem with pure in silico studies is that there is essentially no confirmation towards the in vivo relevance. Drug companies simulate and also test many thousands of compounds in a high-throughput fashion on a daily basis so it is often questionable to what extend individual pure in silico studies without in vivo confirmation do actually add to the bigger picture.
That might be the reason your work was rejected but as said this is a rather general comment and I can neither judge content nor quality of your work so please do not feel offended and read it with a grain of salt.


Scooping and GeneMatcher A colleague has uploaded about a dozen candidate disease genes to GeneMatcher, described here and here and located here,  with excellent success. However last year and again recently he has encountered the case wherein a case was uploaded and a cohort formed, but a few months down the line one of the dropped out collaborators published an n=1 paper in a low tier journal, basically using GeneMatcher to check whether the disease link of the gene was legitimate, to assess the phenotypic range and see what test are worth doing.

Is there documented precedence of this behaviour (preferably citable)? If not,
does anyone have personal experience of this and what is the ratio of incidences over total submissions (I know that is hearsay, but it would be very helpful)?
And also has anyone successfully taken action (say getting user barred)?

I am having trouble finding anything that directly addresses the issue of failed collaborations leading to scooping, but here are some examples of similar behavior:

it is sort of implicit in the discussion here.

This question on academia SE is related but not exactly right.

There are some related stories here.

Here is an article about plagiarism among collaborators, which is quite similar.

I feel like a lot of the discourse around preprinting etc. is similar but again not quite there.

As for punishment, that sounds hard. Scooping is sort of built into science; fighting over precedence is a big deal and has been for a long time. My solution would be to get rid of current academic publication and have everyone's open digital lab notebooks streaming results to publication feeds in real time; that way there's no ambiguity. Preprints and micropublications are sort of a poor man's version of this.
This is not a very complete answer but there may be enough there to help.


Error Code in WGCNA: Dealing with large data sets: block-wise network construction andmodule detection I am comparing single block module detection with 2 block module detection in WGCNA.

Does anyone know what this relates to? I have rechecked and rechecked my codes for both the single and 2 module detection and cannot see an error. Thank you

I get the error code:

Error in plotOrderedColors(dendro$order, colors = colors, rowLabels = rowLabels,  :
  ERROR: length of colors vector not compatible with number of objects in 'order'.

Using the following code:

net= blockwiseModules(Th1_h5, power=6,
                      TOMType = "signed", minModuleSize=30,
                      reassignThreshold =0, mergeCutHeight = 0.25,
                      numericLabels = TRUE, pamRespectsDendro = FALSE,
                      saveTOMs = TRUE,
                      saveTOMFileBase = "PBCreponseTOMB",
                      verbose= 3
mergedColors = labels2colors(net$colors)

plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

bwet = blockwiseModules(Th1_h5, maxBlockSize = 8000,
                        power=6, TOMType = "signed", minModuleSize = 30,
                        reassignThreshold = 0, mergeCutHeight = 0.25,
                        numericLabels = TRUE,
                        saveTOMs = TRUE,
                        saveTOMFileBase = "Th1_NR_R_TOM-blockwise",
                        verbose=3)
bwLabels = matchLabels(bwet$colors,moduleLabels)

bwModuleColors = labels2colors(bwLabels)

plotDendroAndColors(geneTree,
                    cbind(moduleColors, bwModuleColors),
                    c("Single block", "2 blocks"),
                    main = "Single block gene dendrogram and module colours",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang =  0.05)

The problem is that blockwiseModules split your data into (probably 3) blocks, because the default maxBlockSize is 5000 and smaller than the number of genes in your data. If you have enough RAM (16GB is recommended for 12k genes, but 8GB just might be enough), rerun blockswiseModules again with maxBlockSize raised to above the number of your genes. Since you were able to run the function with maxBlockSize=8000, hoepfully you'll be able to do the single block calculation.


Using bioconductor from Python Has anyone used bioconductor from Python? Is there any reason I would choose to use it from R instead of Python? It seems like there is a Python extension for it. Also, is there any reason to use bioconductor over Biopython?

Bioconductor packages are R packages that can only be used directly from python via rpy2 or similar. If you've ever had reason to use rpy2 you'd probably realize quickly that while it works, it's no where near as easy as doing the same thing directly in R.

Biopython is a package for doing a number of useful things, but it can in no way come remotely close to performing even a percent of what Bioconductor packages do. That would be an impossible feat for a single python package.


Extracting all reads from bam file which match read IDs in another file I have a long list of read IDs of interest to me in a file called read_names.txt. it is simply in the format:

m54197_180831_211346/4981510/ccs
m54197_180831_211346/6226723/ccs
m54197_180831_211346/6619607/ccs
...

etc where these are the actual read names from a fastq file. I am then trying to find those in a bam file for which the fastq has been mapped, and I can accomplish this like:

for ID in `cat read_names.txt`
do
samtools view inbam.bam | grep $ID >> read_locs.sam
done

However, this method is obscenely slow because it is rerunning samtools view for every ID iteration (several hours now for 600 read IDs), and I was hoping to do this for several read_names.txt files.

I tried sort of flipping the script a bit and running samtools view first but it only returned the first read ID present in the file and stopped:

samtools view inbam.bam | for ID in `cat read_names.txt`; do grep $ID >> read_locs.txt; done

Am I missing something or is this the only way to accomplish this task?

It is still slow but grep has a -f option to take in a file

samtools view inbam.bam | grep -f read_names.txt > read_locs.txt


How can I subset a reference based on only the first chromosome? I have the GATK FASTA reference for hg19. I would like to subset this based on the first chromosome. Normally what I do is open the FASTA as a text file, and then delete all other chromosomes except chromosome 1.

That's awful.

Picard tools has an option. But this doesn't appear to subset only by chromosome. How could this be done efficiently?

It sounds like you want samtools faidx foo.fa followed by samtools faidx foo.fa chr1 > your_subset_file.fa (or whatever the first chromosome is). The output file is then a regular fasta file subset as you like. You can get any chromosome you want with that. In fact, you can also do regions (e.g., chr1:100-1000, though note that the sequence name in the resulting fasta output will be chr1:100-1000, which may break some downstream tools).


Generate ligands candidates based on protein shape Recent approaches to novel drug design using machine learning (ML) and deep learning, often involve generating hundreds of potential ligands which are later tested by docking with a target protein and recording the resulting binding affinity.

Is it possible to go the other way around by first constraining the ligand generation process to the shape of the protein thus limiting the search space, as opposed to generating ligands of random shapes?

TL;DR: docking is much slower than any ML approach, but the ML approach can be constrained by pharmacophores dictated by the active site.

Side note: Scale

The scale for ligand space exploration is generally is generally several orders of magnitude higher than "hundreds": Zinc DB lists 750 million enumerated compounds, GDB13 has enumerate all possible compounds up to 13 atoms (970 million), Enamine Real is 1.2 billion etc. In an actual experiment, one would use a more restrictive set. As an example of a real project, in Covid19 moonshot project there are 2,800 user suggested follow-ups.

Side note: More than just binding

There are a few requirements for a drug to be effective, and binding is only one of them:

Must follow Lipinksi's rule of five, i.e. it must be membrane permeable and small.
Must bind.
Must be synthesisable. Whereas the experience of someone learnèd is by far the best away to tell, there are several machine leaning approaches that can filter a huge dataset. IBM actually have a fun online tool.
Must not have side effects

So pre-filtering the list of compounds, regardless of docked poses, will happened, be it manually, by cut-offs or by advanced machine learning.

CPU time: docking vs. machine learning

These various steps are generally followed in that order.
Calculating properties (MW, logP etc.) for a given SMILES string to filter based on Lipinski's rule takes a fraction of a second, so that is generally always a first step.

Whether the compound can synthesised or bought is a general left last and done manually. Pharmacology is complicated so is ignored.

Docking a ligand takes time and different approaches have different time scales. A decent algorithm that uses implicit solvent model takes about a minute per core. A short MD simulation takes an 1 hour per core or more. There are docking programs for screening that can be faster but these dock only one conformer of the small molecule against a rigid protein and the results are not of much value. Also setting up docking properly is not easy (cf. steps required for Autodock 4, each with their own caveats). So to recap what needs to be considered that results in the slowness:

Hundreds of conformers of the ligand
Some flexibility of the protein, side-chains repacking only or even backbone changes.
thousands of poses explored.
Optionally. Implicit water is not as good as a TIP3 water model

Machine learning have a much smaller time requirement and is way more attractive for publications. Generally these work by using ligand only information, such as properties, splitting it into pharmacophores, colour etc. This requires a decent dataset, such as empirical results, say from a fragment screen, what compounds bound and which didn't. Often there is a lot of hype with machine learning and the result is something close to something already seen, as well put in this blog post by Pat Walters. This is a side effect of the starting data. The recent Stokes et al. 2020 Cell paper uses AI perfectly because of their staggering automated empirical data collection pipeline.

ML after docking

However, ML can actually be used after docking to improve the quality of the score.
The Boyles, Deane, Morris 2019 Bioinformatics paper investigated ligand-only machine learning to see why do ML methods that rely solely on ligand somehow often work. It is a good paper from a sceptical computational biochemistry viewpoint.
The scorefunction used in docking by the sampler can be classified as forcefield-based, empirical or hybrid. In some cases empirical cases, many terms are ligand based properties (/features) that would normally go into a ML search (say one trained on failed hits and successful hits).
For example Autodock Vina (not 4) has regression derived factors in the final weights of the score

Structural features in machine learning: HotSpots

Information of the protein active site can be used in ranking compounds in a non-docking machine learning approach. Namely, finding out what kind of chemical groups bind preferentially to different parts of the protein's active site, i.e. making a "fragment hotspot map". Specifically, a hotspot map allows the determination of a pharmacophore, which is an abstract generalisation of a group of molecules with both both electronic and steric features —a coarse-grain model, if you will. This approach is also used to detect where the active site actually is in the extreme cases where it is not know (primarily a problem for pipelines, not humans). There are several implementations as described here, the first being FTMap (site, paper), which docks 16 different types of compounds. Waters can play a role and severely complicate matters, hence approaches can be taken to counter that. However, unlike docking protein rigidity is not a big nuisance: this is because a hotspot map is most often used for filtering/weighing by the "color" not the shape of the molecules. But shape based searching is possible.

Shape based scoring

In terms of shape-based searching, this is actually a question of having pre-computed conformers, which can be problematic (but can and is often done on the fly), and shape based scoring —for which there are different metrics, such as the classical RMSD to more complexly weighted scores such as SuCOS.


groupGO parameters explanation (ont and level) I need to use groupGO in clusterProfiler to find functional profile of a list of genes, and I am having trouble finding out what some parameters of the function mean and which I should select for my specific case.

The function is defined as:

groupGO(gene, OrgDb, keyType = "ENTREZID", ont = "CC", level = 2,
  readable = FALSE)

I have mouse genes dataset. What I should select for ont and level parameters? What does ont types: MF, BP, CC mean? What is level? Where I could find this info? It is definitely not in the vignette...

As mentioned in the comments MF is molecular function, BP is biological process and CC is cellular compartment. These are the 3 domains of the ontologies. The level refers to the level in the ontology graph. In the example in that link, pigmentation would be level 1, pigmentation during development level 2, regulation of pigmentation during development is level 3 and so on. You can find more discussion of this in this post on biostars from a year ago.

Regarding the exact settings you should use, I don't personally find CC to be an informative domain. The first couple levels are typically too generic to be useful and they'll include too many genes anyway for you to likely pick up any changes. I'd start around level 3 and see how things go (in other programs you would get all of the levels and domains at once).


How to interpret the GROMACS MD Simulation RMSD Graph? Hello brilliant researchers!
After virtual screening, I performed MD Simulation via GROMACS to check the stability of the protein-top hit complex by following this tutorial http://www.mdtutorials.com/gmx/complex/09_analysis.html

After that, I facing some difficulties in interpreting the resulted RMSD graph as it is quite different from most of the published RMSD graph.

Based on the graph attached, how to interpret it? Is AChE-dexketoprofen complex (Red) more stable (or conserved) than BuChE-dexketoprofen complex (Black)?

The opposite-ish.

Your graph is showing only the ligand, but RMSD is generally used for structures as a whole —so I'll mention about whole protein for future readers.

RMSD is a measure of how much does a structure deviate in position from a reference. Your ligand in the protein in red "wiggles" more than the one in black.

Structure RSMD

If this were RMSD of the protein (and not the index group for the ligand), it would not really mean that the black protein is more stable. Stability is to do with force-field energetics not position. Some protein have two state positions, say MBP is a flexing dumbbell. This would be seen as a jump in the RSMD graph, similar to that seen for the black trace.

However, the command you used was based on gmx rms -s em.tpr -f md_0_10_center.xtc -n index.ndx -tu ns -o rmsd_jz4.xvg from the page you linked. Key here is the index -n index.ndx, which limits it to the ligand. Before this happened gmx trjconv aligned the protein in the trajectory, so that common mistake was not done.

Ligand RSMD

In terms of ligand RMSD, lower deviation/less wobble correlates with more stable binding, but energetics (gmx energy) will tell you how strongly does it bind.

The black trace has a jump, so your ligand in the black trace (which wobbles less) has an alternative binding mode. This is why you are looking at RSMD.
To explore this further you can:

look at the trajectory
in PyMOL using the command split_states twice, once for the first and then the last state, followed by align of the two
If it were the whole protein, looking at the per bond matrix -bm option in Gromacs rms and spotting the change.

3Å??

I strongly suggest checking the trajectory because both have a deviation of 3Å, which for a ligand is not good.

Did the energy minimisation converge properly?
If it is because the proteins themselves wobble substantially, revisit gmx trjconv with a focus on the active site minus the ligand.
Is the ligand binding but leaving one substituent group to wobble?

Energetics

As with docking, the "score" is free energy difference (in kcal/mol). This is not the same as k_off and k_on and their ratio K_d (Dissociation constant). For those you need to do steered MD or have a way longer simulation (not recommended). The reason being is that the ligand has to make its way into the active site. However, you have a secondary binding mode in black, so have two free energy scores, this means that there is a stronger k_on.


correlation coefficient versus DEGs analysis: what's the best approach for low expressed genes? I have 6 folders. Each one contains 7 datasets of a specific type of cancer (RNA-Seq) and 7 datasets of normal tissue (healthy control for that type of cancer). A total of 84 datasets.
I want to investigate my gene of interest in each type of cancer.
I believe the best strategy is to merge the FPKM columns and perform correlation analyses of my gene of interest (GOI) versus all other genes in both cancer and normal tissues. As a result, I will have 14 correlation tables, one for each type of cancer and one for each control tissue. Then, I can investigate the pathways related to the genes most correlated with my GOI in each type of cancer and write my biological interpretation.
Am I in the right direction? Is this a good strategy to answer to my question?
I could also compare the control and cancer samples to investigate the significant DEGs (fold change > 1) and the pathways activated by these DEGs, but my gene of interest doesn't have a strong pattern of gene expression, so I'm afraid it won't be in the list of statistically significant DEGs, that's why I'm preferring the correlation analysis approach.
Any other idea?
The controls were collected from the same patient and processed at the same time. I have the gene counts and can use EdgeR, but what should I do if my GOI doesn't appear among the statistically significant DEGs (fold change > 1)? How can I investigate the role of one single gene in a RNA-Seq dataset if it is not differentially expressed? It can have a role in cancer, even if it is not an important DEG. That's why I thought correlation analyses could give me interesting information.

First of all, RNASeq is extremely sensitive to batch effects.  Are these matched controls processed at the same time as the tumors?  IF they are not, lots of your differences will be caused by batch effect, and not cancer.

I believe the best strategy is to merge the FPKM columns and perform
correlation analyses of my gene of interest (GOI) versus all other
genes in both cancer and normal tissues.

What you propose isn't publishable.  Use the software that people have been using and validating for years.  Use either DESeq2 or EdgeR, or limma if you don't have access to raw gene counts.
If the underlying counts are low, the correlation won't mean much, because your real values could be quite off. If your gene has 2 counts in one sample, and 3 counts in another, that does not mean it really correlates with a gene that has 500 reads in the one sample and 750 in the other. I don't know that you can trick your way around it; you need more reads, or more samples for more power.


Convert CD (cluster of differentiation) into gene symbol I want to convert names according to the "cluster of differentiation" nomenclature into accepted gene symbols. For example, CD3 should map to "CD3a;CD3b;CD3c;CD247" and "CD45RO" should map to "PTPRC". If the surface antigen isn't coded for by a gene, it should map to "" or something similar.

I've seen Obtain CD (cluster of differentiation) marker expression information, but the Abcam table used there does not give the gene symbol -- it gives a list of alternate names. I've seen table 1 of this paper, but it contains errors -- for example, TNFRSF11 should be TNFRSF11A. I could copy this table by hand, or I could download this one, though that's annoying because it doesn't have one row per CD.

Is there a canonical (or just a convenient and reliable) way to do this?

I ended up using option 4 and parsing the columns for "Synonym" and "Previous Symbol".

https://www.genenames.org/cgi-bin/genefamilies/set/471


Finding out-group in unrooted tree Im trying to find the out-group in the alignment my data files file a and  file b both the case i'm getting NA when i run this     pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]

I tried NJ instead of upgma i do still get NA not sure what is the issue ,but at the same time when i used other alignment file with 50 sequence  works fine.

a<-read.alignment("abhi_seq/final_abhi_tempa1.fas", format="fasta")
a.phydat<-as.phyDat(a)
dist.a.phydat<-dist.dna(as.DNAbin(a.phydat))
upgma.a<-upgma(dist.a.phydat)
parsimony(upgma.a,a.phydat)
pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
pars.a$tip.label[pars.a$edge[which.max(pars.a$edge.length),2]]

the fact is that the present upgma tree as a midpoint rooting here, corresponding to node 46.

getRoot(upgma.a)

so for this kind of tree, I recommend you to use the midpoint() rooting function.

pars.a <- optim.parsimony(upgma.a, a.phydat)
pars.a<-acctran(pars.a, a.phydat)
plot(midpoint(pars.a))


Why chose CNN for a variant caller Google released their variant caller DeepVariant which won the highest SNP performance award in the Precision FDA Truth challenge (99.999% accuracy).

From the linked github repo, we see that DeepVariant is a CNN, we provide images of the reads as an input and train the neural network with a stochastic gradient descent algorithm.

To me this is unintuitive. Why use images of reads when you already have the reads data?

It's more "images of reads" rather than actual images of reads. In reality, they're feeding in a pileup, which is just a matrix of letters or numbers (you can visualize this as an image, of course). This is more apparent if you look at one of their example python notebooks, CNNs are popular in the image processing field, but at the end of the day images are just large matrices of values so any time you can convert your input into that sort of format then a CNN can be useful (provided you have sufficient training data and some GPUs to throw at things).


Are codons in RNA layered? Are we misinterpreting RNA codons? I am analyzing nucleotide base-pair patterns in RNA and DNA, and had a thought about RNA and DNA (Let me first state though, I am not a biologist; I am an algorithmatician, employing a sort of recursive sequential neural network).
In any event, I was wondering if we are misinterpreting the way that we look at codons in RNA.  Are we?  Are we missing something when we interpret them the way that we do?  A picture is worth a thousand words, so I made a powerpoint slide of my question and screen captured it.  It highlights the specific question I am asking, which is: are codons layered?
Why is this important to me?  When I look at patterns, I know for a fact that unless it is the end of the sequence (a stop codon) that if I have a U-nucleotide, then the next two nucleotides cannot be AA, AG, nor GA.  Additionally, I know that unless there is a U-nucleotide, the sequence will definitely not end in two or less nucleotides.  This is what got me to wondering if there indeed is layered information in the combination of nucleotides.  Please see the picture.
So if we have (as in the picture) an RNA sequence:
AUGGUCAGUCCAUAA
We might interpret that as 5 codons {Methionine, Valine, Serine, Proline, STOP}
But what if there are actually 13 codons: {Methionine, Tryptophan, Glycine, Valine, Serine, Glutamine, Serine, Valine, Serine, Proline, Histidine, Isoleucine, STOP}

This relates to the concept of the reading frame, which is the point at which you start dividing 3-letter sequences into codons. As you point out, there are multiple ways to divide a string of amino acids into triplets - in fact, there are six unique ways, found by shifting the start letter by one three times, and doing the same but reading in reverse.
Typically, there is one commonly used frame which can be transcribed and turned into protein, called the open reading frame. But it is not unheard of for a single piece of DNA to encode multiple, overlapping genes, depending on where you start reading it. So, the sequence you show cannot encode 13 amino acids in a single string, but it can encode three different strings, one of which has 5 amino acids and the other two which have 4 amino acids each.
See Multiple Reading Frames at Wikipedia, which shows the following example of a single piece of human DNA which encodes two different genes, MT-ATP6 and MT-ATP8.


Seurat Violin Plot: Why do dots align in one row? I am analyzing a rather small scRNAseq dataset (39797 features across 3522 samples) and I have observed that dots in the violin plot are plotted in one line next to each other rather than in a cloud of dots (see plot below):

I am just wonndering what the reason might be for this?
I have followed the Integration and Label Transfer tutorial using SCTransform (https://satijalab.org/seurat/v3.1/integration.html). After the integration process I switched back to the SCT assay, which should contain the normalized counts. For plotting I used the following command:
VlnPlot(..., assay = "SCT", slot = "data"), again specifying to use the normalized counts.
Can someone explain me why this is the case?
Thank you already in advance!

That is because SCTransform (mostly) outputs discrete values, you can check the @slot of your Seurat object.


Tools for quality trimming at 5 prime? I'm looking for a tool that is capable to do quality trimming at 5' end and it is configurable. For example, I can choose the quality trheshold, the read length, etc. Any recommendations?

I'm currently using sickle for quality trimming. However, one of my samples have very low quality at the beginning as well and sickle doesn't seem to be able to deal well with it.

You can have a look at cutadapt. It is capable of quality trimming for both ends as you can read here.


What is the best argument for maximum parsimony method in phylogenetic tree construction? Beside it (maximum parsimony) being compuationally cheap what other good arguments are there for it? Is there any model behind this principle ? Why would one expect this principle to provide right phylogeny in any situation at all ?

I'd say maximum parsimony approaches assume genetic similarity is unlikely to occur by convergence, but likely a result of divergence. If this assumptions holds, finding the tree configuration that minimizes the number of unlikely events makes sense.


What's the difference between VCF spec versions 4.1 and 4.2? What are the key differences between VCF versions 4.1 and 4.2?
 It looks like v4.3 contains a changelog (specs available here) but earlier specifications do not.

This biostar post points out one difference: the introduction of Number=R for fields with one value per allele including REF — can anyone enumerate the other changes between these two versions?

This is easy to check, you can download both specs in .tex format and do diff.

Changes to the v4.2 compared to v4.1:

Information field format: adding source and version as recommended fields.
INFO field can have one value for each possible allele (code R).
For all of the ##INFO, ##FORMAT, ##FILTER, and ##ALT metainformation, extra fields can be included after the default fields.
Alternate base (ALT) can include *: missing due to a upstream deletion.
Quality scores, a sentence removed: High QUAL scores indicate high confidence calls. Although traditionally people use integer phred scores, this field is permitted to be a floating point to enable higher resolution for low confidence calls if desired.
Examples changed a bit.


Is it possible to use SNP heterozygosity as a proxy for Indel heterozygosity? I have estimated genome-wide heterozygosity levels using maximum likelihood and classical substitution model implemented in package atlas. These estimates are way more robust than classical SNP calling and simply counting number of heterogeneous locations because they integrate information from a big window into one estimate. I.e. I am quite sure about heterozygosity levels of SNPs of my genomes.

However, now I would like to know also about heterozygosity of indels and call individual snps using my robust estimates as a prior. GATK allows to specify priors on heterozygosity both of SNPs and indels, which might improve variant calls a lot, because default GATK is parametrised for human (Default priors for heterozygous SNPs is 0.001 and for indel 1.25E-4.). But, I work with insects I know that my SNPs heterozygosity is somewhere between 0.01 and 0.03 (1% - 3%).

So here comes the question, if I know the heterozygosity of SNPs, can I simply guess a reasonable prior for the heterozygosity of indels? GATK defaults for priors are 10x different, should I assume the same relation? Is there a reasonable way how to guess a prior for heterozygosity of indels?

SNPs, insertions, and deletions are generated by different biological processes (see for example here). I would expect the per-generation rates of these mutational processes to vary across species. I would also expect the degree to which each mutational process varies across species to vary as well. I think the answer to your question is "probably not without further information".

However, priors are a feature of Bayesian statistics. In Bayesian inference, it is good practice to redo one's analysis using several priors to establish how sensitive one's results are to the chosen prior. I would suggest that you try to call indels using a range of priors and compare the quality of your call sets to determine how important the choice of prior actually is.


Mapping a list of cells in seurat featureplot I have 209 cells, I clustered them by Seurat to 4 clusters. By Featureplot I am able to track a gene in clusters:

Higher color shows higher expression.

Now, for some genes I want to highlight some cells in Featureplot so that apart from yellow or red colours I want to colour a subsets of cells with another color. I mean I want to map a list of cells in Featureplot or tSNE plot. Let's say I want to know the location of cells 1, 4, 80 and highlight them with another color.
My Seurat object in this link.

Seurat itself beautifully maps the cells in Featureplot for defined genes with a gradient of colours showing the level of expression. Saying I have genes A and B, in excel.

I have coloured cells that express a gene > mean + se, < mean - se or between these values. For instance, for this gene, 36 cells express this gene > mean + se, I want to map these cells in Featureplot or tSNE plot in distinct colour so I can locate them in clusters easily. Something like binary (on off) expression to relative expression.

How I can change the colour of some cells based on my threshold for this gene please?

To color the TSNEPlot, you can generate a new column in metadata with the expression levels (High, low, etc). Then use pt.shape to set a shape for each identity.

To show binary expression based on expression you first have to define the list of cells that are below or over your threshold. Once you have those lists you can use SetIdent() in Seurat to color those groups.

load(file= "~/seuset_16.RData")

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
seuset_h16 <- StashIdent(seuset_h16, save.name = "expr")
seuset_h16 <- SetAllIdent(seuset_h16, id="res.1")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, pt.shape ="expr")
t

high <- c("s1.1","s1.4", "s1.80")
low <- setdiff(colnames(seuset_h16@data), high)
seuset_h16 <- SetIdent(seuset_h16, cells.use = high, ident.use = "high")
seuset_h16 <- SetIdent(seuset_h16, cells.use = low, ident.use = "low")
t <- TSNEPlot(seuset_h16, pt.size = 2, do.return = T, colors.use = c("red", "grey"), do.hover =T, data.hover =c("ident", "res.1))
t


Merging transcriptomes coming from different experiments I'm planning to build a transcriptome by pooling all existing transcriptomes in SRA for a non-model species (which has no reference genome) to study differentially expressed genes and the like. It might worth mentioning that the transcriptomes I want to pool mostly come from different tissues and sampling sites. These transcriptomes would come from five SRA BioProjects plus one project from my laboratoryfor a total of 102 single transcriptomes coming from about 20 individuals (plus three pooled larvae transcriptomes). I'm especially interested in differentially expressed genes in different tissues, on genes specific to different tissues and on the UTRs of differentially expressed genes.

My idea is to build a common reference transcriptome against which I'd re-assemble the individual transcriptomes, and only then I'd get differential expression data etc.

I'm looking for literature (or personal experience) dealing with assemblies of transcriptomes / comparing transcriptomes coming from different datasets in order to understand whether this is a good idea or not and get a grasp of what technical and biological problems I could face. I'm (vaguely) aware of the batch effect and of one way to remove it, but I'd like to read papers dealing with my same problem (different SRA databases, without reference genomes). Up until now I found just projects dealing with microarray data.

As mentioned in the comments, your biological question is going to be vitally important in determining the best way to tackle this issue. However, there are some things that you should generally keep in mind.

Most de novo assemblers will perform kmer normalization. This means that the merged RNA-Seq data will be normalized across all libraries. This is almost certainly not what you want because you'll be randomly selecting reads from all libraries without regard to their treatments. So, if you're going to allow for kmer normalization, you should normalize prior to merging the libraries so that you have equal representation of the reads by library. If you don't normalize, you're going to have a very fragmented assembly, and it will take a long time to generate.
De novo assemblers using short reads assume all identified exons with junction support may be used to generate splice variants with equal probability. I'm not entirely sure how to best phrase that, but here's an example. Suppose in a cold-treated library that Gene A generates a transcript using exons 1-2-3, but a heat-treated library generates transcripts using 3-4-5 and 2-3-5. It may be that merging the libraries results in transcripts generated containing various permutations of those exons (e.g. 1-2-3-5, 2-3-4, etc.) that aren't generated in any of the biological settings but can no longer be distinguished. This will lead to false positive transcripts which will detrimentally affect expression quantification.
Merging assemblies is difficult. So, to be true to each treatment, it makes more sense to generate individual assemblies (which will also result in library-specific kmer normalization). StringTie essentially does this using reference genomes. You assemble using each library (all bio-reps for a treatment), and then merge after the fact and re-quantify. This is the best approach. However, it can be difficult and computationally taxing. Some programs like CD-HIT (specifically cd-hit-est) are typically used in these cases. Keep in mind that CD-HIT will drop sequences with high similarity, not merge them. This may lead to loss of some uniquely assembled sequence. Other approaches involve using programs like CAP3 (a DNA assembler) to treat the transcripts like longer reads for assembly. Depending on the parameter choice, this might work, but you will risk generating chimeric transcripts.

How you move forward with the assembly process will again depend on your biological questions. I would tend to promote that you tackle the questions one at a time rather than using a single transcriptome to address them all, but there may be circumstances in which merging assemblies makes more sense.

In your specific case, I would perform multiple assemblies using a variety of the approaches described above and then annotate the assemblies (ex. using Trinotate) to see how their quality metrics compare. I would discourage generating a single transcriptome and then re-assembling individual transcriptomes from that assembly -- the single transcriptome may have misassemblies that would then propagate to the individual assemblies, and this process would generally be messy/difficult.


Machine learning using protein-sequences I'm participating in a bioinformatics machine-learning seminar at my university. The main task is predicting binary classification of protein-protein interactions using sequence data as input.

One of the subtasks is familiarization with the dataset and presenting the dataset. Now I'm wondering which additional information I can get out of the sequence data.

I want to start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact. The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”. Species are completely mixed.

Is there something interesting that I can additionally include in my analysis? I'm used to exploratory data analysis though this is different due to the nature of the protein sequences.

There are two types of features which can be extracted from protein sequences.  It is not possible to known which, if any, of these features would be useful in classification tasks.  It may not be possible to build a classifier at all, or it may be very straight forward.  To know this, a feature selection technique must be used, such as e.g. forward feature selection.

Sequence based features

Sequence based features are derived from sequence-only information.

The most common used features are the PseAAC (pseudo amino acid composition) style features, including:

Average Position of each amino acid - 20 features
Average Distance between each amino acid - 20 features
Motifs - 20 features, 20*20 features, and/or 20*20*20 features (depending on motif length)
OAAC (Overall Amino Acid Composition) - 20 features
SAAC (Split Amino Acid Composition) - 60+ features
AAIndex Composition - 1 to 567+ features
PCC (Physicochemical Composition) - ~13 features
PSSM encoded as 1, 20, 200, 400 or 600 features (i.e. Psi-Blast output using NR or swissprot)
Average predicted residue burial/exposure level
Average predicted sequence entropy

Structure based features (predicted)

Using the sequence, you can either predict the 3d structure using homology, ab initio, or directly use a homologous structure.  You can then extract predicted 3d features, such as:

3d size (i.e. max of x,y,z)
Calculated ASA/SASA
3d contacts (type, frequency, etc)
any other local/global environmental information
cath/scop classification
dssp/stride composition


Understanding some of the computational bottlenecks of Covid-19 research I am a researcher in high-performance computing (with very little bioinformatics background), and I am trying to understand what are the current biggest computational bottlenecks of software used for research on the topic of Covid-19 (testing, vaccine development, drug identification, etc...).

I know this is a very broad question, but I am looking for the biggest common denominators. Said differently, if you could pick a few operations / algorithms / data structures to speed up Covid-19 research and they would magically run 1000 times faster, what would they be?

Note: links to books, papers, and scientific literature are very welcomed.

I can only speak of drug design (and even then I am terrible at turning down the jargon).

In the case of drug design, this is pretty much plan C. Namely, none of compounds that entered clinical trial at the start of the year work (let's call this plan A) and none of the vaccines that are entering now clinical trial work (let's call this plan B although vaccination prevents so is a better solution). Clinical trials take a lot time as they need to be thorough (cf. the infamous thalidomide from the days with more lax testing)

There are several bottlenecks in drug design. These are technical based on our imperfect understanding of things, which results in less than ideal models. CPU time is not really an issue. Empirical validation is required at many steps along the way. This is greatly increased by automation. The XChem team at Diamond light source have done a fragment-screening against the SARS-CoV-2 main protease at an insane pace: things are changing fast!

Drug screen

There are many ways to make a drug. In one way, a test system was set up (a reporter assay) that gives you a measurable read-out, say mm of clearance of bacteria on a Petri dish (viability) around a cellulose dish soaked in different separated natural extracts —and you might find that Penicillium chrysogenum has a compound of interest. What you are doing is drug screen. Screening a large panel of compounds (of varying purity) does not necessarily require knowing the protein target.

Rational drug design

Fast forward to this decade and it becomes possible to predict the free energy release upon binding of a small molecule against a protein —mostly.

This requires a solved 3D structure of the protein by X-ray crystallography. Crystallising a protein is a tricky business, as a pure protein sample has to come out of solution in a folded form to form a neat lattice as opposed to just making unfolded gunk. What salt conditions work? Anyone's guess. Some have done deep learning on it, but the best solution is to have a robot try a very large range of combinations.

Once you have a protein structure you can predict what fits into the part that you want to target, say an active site or an interface. This is not done by simple cartesian geometric shape fitting like LEGO bricks as there are many molecular forces involved. So a mathematical model of these forces is made, called a force field. Many have been made, a lot of study has gone into these, but they are not perfect. It is all made worse by the fact that we cannot really model water properly —most models use a proxy (implicit solvent) or use non-polarisable water.
This is all while made worse by the fact that a protein is not static but wobbles (following the same mathematical models). So a ligand does not slot in really like a key in a lock, but has an induced fit. Modelling all this is accurately is problematic, but one can still try: calculating where exactly a ligand binds to a protein is called for which one gets energy difference upon binding. And one can try an inordinate amount of small molecules —Zinc DB contains 230 million purchasable compounds. This is called a virtual screen. One can even start with the natural substrate and "weaponise it" by adding chemical groups that do a chemical reaction that breaks the protein. However, a metabolite is very different from a drug (cf. Lipinski's rule of five) as the latter must cross membranes etc. But generally the accuracy is not great without multiple steps of empirical validation.

To improve the results from a virtual screen the most common approach is to do a fragment screen to guide it. This is to automatedly soak protein crystals with a large library of small compounds and see which crystal contains a bound hit. The hits can be expanded into bigger fragments in a virtual screens and the short list of compounds tested (these are not cheap) by soaking and crystallising again.
In the case of Covid19 virus, the cool thing is that this whole process is open in the aforementioned current campaign (Diamond is a government owned synchrotron and not a big pharma). What is really nice about this is that for the second cycle comp-chemists from everywhere have done their own guided virtual screen and have submit compounds (which you can see publicly).

It should be said that drug screens that are not rational design are still a thing. The recent antibiotic-discovering deep-learning paper in Cell]6 is such an example, although it should be said that it was a feat not only in AI but also in automated empirical testing.

Future

Automation via better and open software will definitely improve. With a small drive possibly from the dissemination of easy to use and to code off-the-shelf electronics, such as Rasperry Pis or Arduinos. But mainly with more complex machines such as the Opentron system (a python-controlled robot).

But indubitably where the greatest improvement can come is better mathematical models. Possibly even with deep learning. After all researches were amused last year-ish (2018 competition) when DeepMind "won" the protein 3D structure (fold) modelling competition (CASP13) and there have definitely been some really nice results in the field of docking with hybrid models. So deep-learning may still surprise us yet...

Lastly, open science will be a strong driver: using premade packages or tools on GitHub and code snippets from StackOverflow etc. will accelerate these.


Method to detect genome doubling CNV is usually detected using either read depth or SNP VAF. However, neither of these methods can detect CNV in case of genome doubling (assuming genome doubling is a kind of genome scale CNV comparing to normal tissue which is diplod). We thought of other platform such as SNP array and we did not think they can detect genome doubling as well. Can any one share some comments on the method for such event? Thanks in advance.

In a tissue or cell line that used to be diploid and has undergone genome wide doubling, you expect more variants having observed allele frequencies that would also occur in tetraploid organism. That is, observed VAF distribution will have (possibly small) peaks at 0.25 and 0.75, in addition to the diploid peaks at 0, 0.5 and 1. Compare the distributions in the tissue or cell line of interest to those of the control, that is, known diploid. Also if it is a doubled male genome, the sex chromosomes will show a new peak at AF 0.5, that is, heterozygous X or Y variants. This is just an idea, I have never tried to look for such events.


Least present k-mers in the human genome What are the least present k-mers in the human genome at different sizes?

Starting with k=4 and going up in size until k=10, what are the k-mers least seen (or not at all) in the human genome? I am only interested in the reference human genome, so I am discounting SNPs/Indels in the population.

If this is not pre-calculated somewhere, what tool is recommended to use that starts with the GRCh38 reference as input?

You can use the Jellyfish software to calculate the k-mer profiles up to length 31.

From the instructions in the user guide:

The basic command to count all k-mers is as follows:

jellyfish count -m 21 -s 100M -t 10 -C reads.fasta

To compute the histogram of the k-mer occurrences, use the histo subcommand (see section 3.1):

jellyfish histo mer_counts.jf

To query the counts of a particular k-mer, use the query subcommand (see section 3.3):

jellyfish query mer_counts.jf AACGTTG

To output all the counts for all the k-mers in the le, use the dump subcommand (see section 3.2):

jellyfish dump mer_counts.jf > mer_counts_dumps.fa


Using sickle for quality trimming I am using sickle package from bioconda to trim my reads from the Illumina Miseq.

When I tried following code, it failed to run

sickle pe -t illumina -f file1   -r file2  -o output1  -p output2  -s combined_output -q 27 -l 25

When I replaced illumina with sanger, it works. But my reads are from the illumina. Could anyone clarify if using this option is acceptable ? Also, what is the optimum value of -q and -l

It'd be rare to run into Illumina encoded quality scores any more, I doubt there are machines running that still produce those. Everything produced at the moment uses Sanger encoded (phred + 33) scores.

There is no optimum value for -q and -l, though reasonable values would typically be <=5 for -q and ~5 for -l.


NCBI - edirect suite to download all genome sequences associated with a query - troubleshooting Not sure if this is allowed but I cannot think of a better place to ask. I am attempting to download all genomic sequences from refseq associated with the query Peptostreptococcaceae. I have been been using the following command (based on some reading of the manual):

esearch -db nuccore -query "Peptostreptococcaceae" | efilter -molecule genomic -source refseq | efetch -format fasta

But I keep getting the following error

Unknown option: molecule
Unknown option: source
Must supply -query or -days or -mindate and -maxdate arguments on command
line Db value not found in fetch input No do_post output returned from
'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=nuccore&term=Peptostreptococcaceae&retmax=0&usehistory=y&edirect=3.80&tool=edirect&email=lgunn@sonic15' Result of do_post http request is .............
(a lot of text which I assume is .html or .xml related)

The issue here was I was using an older version of eDirect. There is no command to check version or updates. Your best bet is if you experience problems uninstall and reinstall.


How to map selected genes to Metabolic pathway Maps I have a selected Arabidopsis Genome Initiative (AGI) list for RNA seq and proteomics data, how can I map them to metabolic pathway maps to vilualize(e.g. TCA cycle / FA / Photosynthesis) in KEGG or AraCyc or any freely available pathway maps?

Could someone share an example if possible?

KEGG has a webtool for this, KEGG Mapper.

You can just copy-paste your gene ids in there and see how your pathways are colored.


How to validate that BAMs have been downloaded correctly? I currently have several hundred BAM files which were downloaded by someone else. These have remained untouched---before working with them, I would like to double-check that these BAMs have been fully downloaded.

I don't MD5 Checksums to look at.

In order to validate the BAMs, I would use something like Picard's ValidateSamFile, https://broadinstitute.github.io/picard/command-line-overview.html#ValidateSamFile

java -jar picard.jar ValidateSamFile \
      I=input.bam \
      MODE=SUMMARY

I could write a bash script, and run through these BAMs, and then read through the Picard output. That is checking whether the BAM is valid though, not whether it was downloaded fully (e.g. the BAM could be downloaded correctly, but be inherently flawed in some other way).

If the BAM is downloaded fully, isn't there some "end of file" bytes I could quickly check? Are there other more efficient methods?

samtools quickcheck is all you need. From the manual:

Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least one target sequence and then seeks to the end of the file and checks that an end-of-file (EOF) is present and intact (BAM only).
Data in the middle of the file is not read since that would be much more time consuming, so please note that this command will not detect internal corruption, but is useful for testing that files are not truncated before performing more intensive tasks on them.
This command will exit with a non-zero exit code if any input files don't have a valid header or are missing an EOF block. Otherwise it will exit successfully (with a zero exit code).