How can I extract gene names for a metabolic pathway from KEGG?

I whipped this little script up using the KEGG API:

#!/usr/bin/env python3
import urllib.request
import re
import sys
pathway = 'hsa00010' # glycolysis
url = "http://rest.kegg.jp/get/" + pathway
with urllib.request.urlopen(url) as f:
    lines = f.read().decode('utf-8').splitlines()
    want = 0
    for line in lines:
        fields = line.split()
        ## The list of genes starts here
        if fields[0] == 'GENE':
            want = 1
            ## The line with GENE is different
            print(fields[2].rstrip(';'))
        ## We reached the next section of the file
        elif want == 1 and re.match('^\S', line):
            sys.exit();
        ## We're still in the list of genes
        if want == 1 and len(fields)>1:
            print(fields[1].rstrip(';'))


Performing t.test to obtain p.values for each gene

You can simply apply() t.test() to your matrix. In general, though, I expect lmFit() from the limma package will be simpler to use.

Here is an example with t.test():

m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data
apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value)

Neither t.test() nor lmFit() should be expected to produce the same results as excel, since excel is inferior in every way and should never be relied upon for any statistical analysis (plus lmFit is doing something much much more complicated than a simple t-test).

As an aside, FPKMs are highly problematic for use in statistics. If possible, I strongly encourage you to get either per-gene counts or per-gene estimated counts and to use those in limma or another differential expression package instead, as the results will be much more reliable. I presume, though, that you're using TCGA data, so you may be stuck with FPKMs.


getting KeyError on certain residues in BioPython: ATOM vs. HETATM lines

You can find the answer in The Biopython Structural Bioinformatics FAQ:

What is a residue id?
This is a bit more complicated, due to the clumsy PDB format. A
residue id is a tuple with three elements:

The hetero-flag: this is 'H_' plus the name of the hetero-residue (e.g. 'H_GLC' in the case of a glucose molecule), or 'W' in the case
of a water molecule.
The sequence identifier in the chain, e.g. 100
The insertion code, e.g. 'A'. The insertion code is sometimes
used to preserve a certain desirable residue numbering scheme. A Ser 80
insertion mutant (inserted e.g. between a Thr 80 and an Asn 81
residue) could e.g. have sequence identifiers and insertion codes as
follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering
scheme stays in tune with that of the wild type structure.

The id of the above glucose residue would thus be ('H_GLC', 100, 'A').

So in your case it would be ('H_MSE', 85, ' ').


The definitive answer is the sequencing institute/lab. They know what protocol/chemistry etc. was used.
If you don't have access to that a number of tools check for known adapter sequences. Run e.g. FASTQC, which will tell you the proportions of adapter sequence found. Tools like TrimGalore can also autodetect the most common adapters.
UPDATE:
I checked 2M reads from one run of your dataset and it seems to me that it was uploaded already preprocessed:

Check all your samples, but I don't see a need to do any preprocessing if all of them look like this.


If you have your counts in a data.frame called counts, something like this might work:

filtered.counts <- counts[rowSums(counts==0)<10, ]

For example lets assume the following data frame.

> A <- c(0,0,0,0,0)
> B <- c(0,1,0,0,0)
> C <- c(0,2,0,2,0)
> D <- c(0,5,1,1,2)
>
> counts <- data.frame(A=A, B=B, C=C, D=D)
> counts
  A B C D
1 0 0 0 0
2 0 1 2 5
3 0 0 0 1
4 0 0 2 1
5 0 0 0 2

And now you would want to remove the genes with at least 3 zeroes or more:

> filtered.counts <- counts[rowSums(counts==0)<3, ]
>
> filtered.counts
  A B C D
2 0 1 2 5
4 0 0 2 1

Every gene with 3 or more zeroes is now removed.


How can I specify an email address in the SLURM snakemake json configuration file?

As noted by the comments and the SLURM manual (not so much bioinformatics related?)

You need the following terms in your batch/queuing command:

--mail-type=ALL
--mail-user=you@email.com

I used "ALL" which will send you all reports. If you just want the finish message then you would use "END" and "FAIL" like such:

--mail-type=END,FAIL

Here are the specifics from the manual if you wish to make a different selection:

--mail-type=

  Notify user by email when certain event types occur. Valid type values
  are NONE, BEGIN, END, FAIL, REQUEUE, ALL (equivalent to BEGIN, END,
  FAIL, REQUEUE, and STAGE_OUT), STAGE_OUT (burst buffer stage out and
  teardown completed), TIME_LIMIT, TIME_LIMIT_90 (reached 90 percent of
  time limit), TIME_LIMIT_80 (reached 80 percent of time limit),
  TIME_LIMIT_50 (reached 50 percent of time limit) and ARRAY_TASKS (send
  emails for each array task). Multiple type values may be specified in
  a comma separated list. The user to be notified is indicated with
  --mail-user. Unless the ARRAY_TASKS option is specified, mail notifications on job BEGIN, END and FAIL apply to a job array as a
  whole rather than generating individual email messages for each task
  in the job array.

--mail-user=

  User to receive email notification of state changes as defined by
  --mail-type. The default value is the submitting user.


Given a transcription factor, what genes does it regulate?

MSigDB has a collection (C3:TFT) of gene sets corresponding to transcription factor targets.

Harmonizome has functional terms for genes extracted from over a hundred publicly available resources.


Download or create VCF from Human Genome Diversity Project (HGDP)

I presume you are looking for the recently sequenced WGS data rather than the genotype data? If so, the different files can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/
Each chromsome has a seperate vcf file with a corresponding .tbi index file.
The gvcfs can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/gVCFs/
These files are more compressed and smaller in size than the regular vcfs.


You are right that a single molecule in a single position is either methylated or not methylated. However:
First, assuming your organism of interest is diploid (or of higher ploidy) one of the chromosomes could be methylated, the other not. That would give you a level of 0.5 and can be found in imprinted regions (where the paternally inherited chromosome is methylated and maternal not, or vice versa).
Second, except if you use a pure cell population of a single cell type it is likely that there will be cellular heterogeneity in your sample, e.g. from different types of blood cells or neurons vs glial cells in the brain. Some of these may have methylated CpGs in this locus (depending on cell-type-specific gene activity) while others can be unmethylated. Depending on the relative proportions of these cell types the methylation level can be anything between 0 and 1 for a mixture of cell types.
Thirdly, sometimes CpGs get grouped together. For example, 5 nearby CpGs are combined to give a single methylation level for this locus, and then the level can be 0; 0.2; 0.4; 0.6; 0.8 and 1.
Lastly, there could be technical artifacts or technical reasons why a position was not called as methylated (insufficient quality for a single read?), leading to deviation from 0 or 1.


If you're in the RepeatScout directory you need to execute the Perl script with a ./ prefix, like so.

cat od_1000_repeats.fa | ./filter-stage-1.prl > od_1000_repeats_filt1.fa

If you're not in the RepeatScout directory, you need to provide the full relative or absolute path when executing the script. Alternatively, you can put the RepeatScout directory in your $PATH environmental variable.

These are all standard UNIX/Linux concepts that many bioinformaticians take for granted, for better and for worse. If you're lost in all the terminology, a crash course in UNIX basics is probably where you would want to start.


How to interpret PCA output statistically and biologically?

All that plots like this are telling you is that there are some genes that contribute more to the variability seen between your various samples than others. In an ideal world these genes will also be differentially expressed between your groups, but since we don't live in an ideal world that might not be the case. In general, the longer the line and the closer it is to the X axis the more it contributes to the overall variation (since it's contributing more the PC1, which is 53.5% of the variation in your data).

Since you've asked in the comments for a brief "results and description" section style summary of such a PCA I'll include something brief below. Note that one should not over-interpret PCA plots. The axes are (mathematically) orthogonal and don't necessarily correspond to anything biologically relevant.

  We performed a PCA on the variance-stabilized counts to check for batch effects and overall clustering of the data. As can be seen, the "3T" and "5T" groups cluster together along the first principal component, while the "0T" and "1T" samples cluster on the opposite side. This suggests that these groups are more similar to each other, as further evidenced by ...

Assuming you had used the vst() function in DESeq2 (I note you've tagged your post with RNAseq, so perhaps you're using DESeq2) then something like that above would be a reasonable description.


finding unique headers in a fasta file using linux command line

If you just want the number of unique headers, you can do this:

grep '>' reference.fasta | sort | uniq | wc -l

If you want a list of the unique headers, you can do this:

grep '>' reference.fasta | sort | uniq

If you want a histogram of how many times each header occurs, you can do this:

grep '>' reference.fasta | sort | uniq -c | awk '{printf("%s\t%s\n", $1, $2)}'


Standard Cutoff for Moderated T-statistics

You're misinterpreting the moderated T-statistic, it's basically the fold-change divided by its variance. The p-value comes directly from that, so if you filter by moderated fold-change you're just setting an unknown (unless you go through the trouble of figuring it out) p-value threshold. Instead, do exactly as you've been doing and set appropriate (adjusted) p-value and fold-change thresholds.


BED format can refer to the ".txt" file format that stores chromosomes and positions of genomic features as well as some other things like name. bed format can also refer to the plink format, which belongs to bed, bim, fam trio. The plink bed format stores all the genetic variants, but the program needs to refer to bim file for their positions.
Genome liftover obviously only needs the coordination of each SNPs, and it does not need individual genotypes at all. In fact, the bed format in PLINK does not provide any useful information. Knowing this means that providing the plink bed file is wrong. And this is pretty clear from the error message.
Moreover, since the three plink file always needs to be used together, most programs refer to them as plink file format. So whenever people say BED format, it is the txt file format.


data = pd.read_excel ("2019-02-27_161601 AtWAKL8 different version expressions.xls", sheet_name='Results').fillna(0)

data.to_csv('df1' + '.csv', index=True)

df1 = pd.read_csv ("df1.csv")

This doesn't make much sense ... you read in Excel to a dataframe, you then read out to csv - the dataframe doesn't have an index, but you include an index and then read it in again. This could do more harm than good.

Anyway your supervisor is right it is very easy to parse this data using pandas, although there will be more expertise doing this via R.

I don't want to peer through the gloom of a CSV output, but I would like to highlight how easy it is to parse this stuff into a dataframe. The dataframe 'data' is fine.

Parse You can simply remove the top rows of bumf data = data.iloc[6:] (I might have miss counted, plus see below),
Headers You will probably have to remove the "headers" row to a list and
then introduce back onto the dataframe using data.columns =
['header1', 'header2' ...] ... rather this is the list (array in Perl) which you can automatically generate using a magic to_list type pandas command which I forget, its not list - its cuter
Columns You next might use drop to remove the columns you don't want (that might be alot of columns) or else make a subset of
the dataframe into a new dataframe, you can also use iloc for
columns too.
Index You then assign A1, A2, ... etc  as an index data.index = ['A1',
'A2', ]...etc..  you can reshape column 3 to automate this
to_numeric You might have to assign the CT values as floats
    or better to_numeric
Extrapolation You would then define a subroutine (or whatever they call them in
Python) so perform the extrapolation and use the apply command on
your dataframe and assign it to a new column, where the input is your CT values.
Output Finally, issue your to_csv command to_pickle is also a good idea because you have a df that is worth importing at a future date.
Plotting this is easy using data.plot which appears to be some sort of inbuilt matplotlib, keep in mind seaborn DONT try using R's ggplot2 within Python using rpy2 .... its a NIGHTMARE to get working and doesn't support the latest ggplot2 (!!!). R
s ggplot2 is stronger  than seaborn or matplotlib (providing you are working within R), but these two are good for quick graphs.

If I had time I'd write the code for you, but I'm afraid I'll have to leave it there. I'm still in debt with a Perl efetch script I promised a previous OP.

In summary, Python pandas is fabulous stuff for this purpose, I love it. Doing this via Perl and "subprocessing" inside R... its doable, well in this particular case I suspect its just as easy, but generally its alot more code and often fiddly.


bedtool intersect to compare .vcf file to reference genome

If you want to filter a VCF file, you might want to give VCFtools a try.

Here is the documentation : VCFtools Manual

I am not a specialist of this tool but you could do as the following. First build a bedfile with the chromosome names and positions of your genes (start and stop) such as this :

chr   start   stop   (name)   (score)   (strand)

Then it should work with this command line :

vcftools --vcf file.vcf --bed positions.bed --out filtered.vcf

The output will be a VCF filtered on the provided regions


Alignment with arbitrary number of mismatches or gaps

This is not something mappers will be good at, because you're going to hit most of the genome. Depending on how many reads you have it may be better to work in reverse, iterating over all the kmers in the genome and working out whether or not they are fewer than 7 mismatches (plus indels) away from any reads:

for gp in genomePositions
    genomeKmer = genome[gp:(gp+23-1)]
    for read in reads
        errors = 0
        for rp in 1..23
            if (genomeKmer[rp] != read[rp])
                errors += 1
        if (errors < 7)
            print gp


Can we describe synthetic biological sequence data as belonging to a "species" or "taxon"?

Phylogenies can be made with/of non-species groups. I think "taxon" is perfectly fine. You might also consider the use of the "leaf / leaves" or "terminal" if you want to emphasise the tip nature of the organisms. Even "sequence" would be fine.

Ninja edit to my own suggestion: since you're making a phylogeny from sequences of unknown context or relationship (i.e. you can't say that they're sequences from a "species" or "OTU" or from anywhere), agnostic terms like "taxon" or "sequence" are more correct.


How are UMIs used to dedupulicate in Drop-seq tools?

UMI and barcode are the same but read2 is different, would they be
  collapsed to 1 UMI count or 2?

If the two reads share cell barcodes and UMI sequences and align to the same gene, (even if they are different sequences) then they originate from the same RNA molecule, and should only count as 1.

https://support.10xgenomics.com/single-cell-gene-expression/library-prep/doc/user-guide-chromium-single-cell-3-reagent-kits-user-guide-v3-chemistry


What datasets are available out there for prediction based on DNA sequences?

Good question and the Janggu link is cool. Thanks!
Kaggle is/was heavily image recognition. There isn't an MNIST type data set in genomics, in case anyone is wondering this is a classic natural language processing (NLP) type analysis to identify handwritten single digit numbers.
There has been several FDA precision challenges, one being described here:
FDA precision challenges
This could provide something towards a 'standard'.
A key issue remains the accurate vectorisation of sequence and whether you think that is/was satisfactory in my personal opinion.
There is an old review article here


How to correlate two zero inflated bedgraph-like signals?

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

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


What you are looking for is "Differential accessibility analysis", you are using ATAC-seq, hence "accessibility". You will need to "call" and then "normalize" "count"(*) peaks, after these, you can use DESeq2 or edgeR for the "differential" part.

Here is a nice Nextflow workflow that addresses all of the steps of a differential accessibility analysis of ATAC-Seq data.

* featureCounts can be used the identify peaks, see this and this.


Make a fasta out of mapped reads without taking into account the reference

vSNP outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  vsnp_merge_vcf_into_fasta.py.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".

vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

Use *zc.vcf output to merge into FASTA

vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf


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.


How to tell which channel is the reference channel in two channel (red green) array data?

No, you can't get this from the data itself. If there is a control/reference sample on one channel, you need someone to annotate which channel it's on.


Network comparison of single cells (from sequencing data)

You can't build a network of a single cell only with the expression of a single cell. You either need previous known interactions or pathways or you need to use several cells/samples.
If you use previous known information, you can use pathway information, otherwise you can group some cells and use something along the lines of WGCNA to find a scale-free network.

In the first case one can add the expression upon it and see what is happening, what edges are kept and which ones not...

Once you have the network for each cell you can compare the metrics for each cell. Some metrics to look at are the centrality of each gene and the betweneess.
In a more general way, you can look if it follows a scale-free topology or other type of network structures, or how many nodes are present, how many edges are present...


Integrative analysis of omics studies using machine learning

The steps you describe are correct. For step 2 it is usually normalized to mean 0 and variance 1. However the "machine learning" part is important.

Having several samples being technical replicates will make the integration task easier. However, you have too few samples to make any good prediction. At most I would describe it as an exploratory analysis.

I recommend using regularized canonical correlation analysis for this kind of multi-omics analysis. The method seeks which are the features of one dataset that correlate (or covariate, depending on a parameter) with the variables of other datasets. Then you could use those features to predict the position of the other.

There are several implementations of this methods, the original in the RGCCA package, and mixOmics for other procedures with cross-validations also STATegRa has other related methods.


SnpEff gene annotation with @ such as PCDHB@

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

It's also mentioned in the GeneCards documentation:


what is the difference between local identity and global identity for homolog finding?

If you want a purely alignment based approach to begin to infer homology (note I say infer, as we stated in the comments, this isn't sufficient to say 2 sequences are homologs), I would use a Multiple Sequence Alignment (MSA) instead of a pairwise like you've done so far, and include as many similar proteins as you can from BLAST.

The choice of algorithm for the MSA probably won't matter too much, as you will be looking for conserved regions anyway.

Once you have your MSA, you can easily spot conserved domains. If anything is known about this protein structurally, it's quite trivial to work out whether a conserved region is a functional domain (active site, that kind of thing).

With this information, you would be somewhat justified in calling these proteins structural homologs at least if they share domains, but, to reiterate, without doing phylogenetic analyses, you cannot call them true homologs.


Group genes by functional categories suming expression values

As suggested in the comments, the pandas module is quite convenient for this type of work.

Here is how you could do (assuming you have tab-separated values as input, and you want similar file format for the output):

import pandas as pd

counts = pd.read_table("counts.tsv")
# At this point, counts looks as follows:
#       gene  s1  s2  s3  s4  s5 Category
# 0   name01   0   2   2   0   0        A
# 1   name02   3   0   1   0   0        A
# 2   name03   0   2   1   0   0        A
# 3   name04   0   0   1   0   0        B
# 4   name05   5   0   1   0   0        C
# 5   name06   1   0   0   0   0        D
# 6   name07   2   0   0   0   0        D
# 7   name08   0   0   3   0   0        E
# 8   name09   1   0   0   0   0        F
# 9   name10   1   0   0   0   0        F
# 10  name11   3   0   0   0   0        F

summed_by_cat = counts.groupby("Category").sum()
# At this point, summed_by_cats looks as follows:
#           s1  s2  s3  s4  s5
# Category
# A          3   4   4   0   0
# B          0   0   1   0   0
# C          5   0   1   0   0
# D          3   0   0   0   0
# E          0   0   3   0   0
# F          5   0   0   0   0

# Change the index name to get the desired first column header in the output:
summed_by_cat.index.name = "Cat"

# Write to a file
summed_by_cat.to_csv("by_cat.tsv", sep="\t")

And resulting the file contains the following:

Cat s1  s2  s3  s4  s5
A   3   4   4   0   0
B   0   0   1   0   0
C   5   0   1   0   0
D   3   0   0   0   0
E   0   0   3   0   0
F   5   0   0   0   0

You could also load your table using counts = pd.read_table("data.tsv", index_col="gene") and counts would have the gene names as index instead of default numbers:

        s1  s2  s3  s4  s5 Category
gene
name01   0   2   2   0   0        A
name02   3   0   1   0   0        A
name03   0   2   1   0   0        A
name04   0   0   1   0   0        B
name05   5   0   1   0   0        C
name06   1   0   0   0   0        D
name07   2   0   0   0   0        D
name08   0   0   3   0   0        E
name09   1   0   0   0   0        F
name10   1   0   0   0   0        F
name11   3   0   0   0   0        F

This doesn't seem to affect the final outcome.


Using pysam with cython: htslib/kstring.h not found

Based on a suggestion by Devon Ryan, I searched how to set C_INCLUDE_PATH for cython, and found this issue. It refers to a (currently not working) feature present in cython's documentation, that should let libraries inform cython of the correct localisation for the header files in the module's setup:

ext_modules = cythonize(
    "bam25prime/libbam25prime.pyx",
    include_path=pysam_get_include())

As noted in the issue discussion, this seems to have no effect, but a workaround exists, consisting in using the Extension class to store the header localization information.

So the general answer is: Some python libraries (including pysam) provide a function that returns the paths to their include directories. This can be used to inform cython, but has to be done "manually".

My setup.py now looks as follows:

from setuptools import setup, find_packages
# If you have .pyx things to cythonize
from distutils.extension import Extension
from Cython.Build import cythonize
from pysam import get_include as pysam_get_include

name = "bam25prime"
__version__ = "0.1"

# https://github.com/cython/cython/blob/master/docs/src/reference/compilation.rst#configuring-the-c-build
extensions = [
    Extension(
        "bam25prime.libbam25prime", ["bam25prime/libbam25prime.pyx"],
        include_dirs=pysam_get_include()),
    ]

setup(
    name=name,
    version=__version__,
    # [description, author, email, licence, ...]
    packages=find_packages(),
    scripts=["bin/%s" % name],
    ext_modules = cythonize(extensions),
    # ext_modules = cythonize(
    #     "bam25prime/libbam25prime.pyx",
    #     include_path=pysam_get_include()),
    install_requires=["pysam", "pybedtools"])

Then, python3.6 setup.py build_ext succeeds.

And for whoever might be interested, here are the other files of the project:

bin/bam25prime is a link to bam25prime/bam25prime.py, which looks like this:

#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
"""
When executed, this module collapses reads in a bam file to their 5-prime end
(reducing them to length 1) and writes the result in bed format.
"""

import argparse
from pybedtools import BedTool, Interval
from pysam import AlignmentFile
from .libbam25prime import collapse_ali

# pybedtools.Interval | pysam.AlignedSegment
# chrom               | reference_name
# start               | reference_start
# end                 | reference_end
# bool(int(strand))   | is_reverse
#
# Interval(reference_name, reference_start, reference_end, strand=strand)

def collapse_and_sort(alis):
    """
    Collapse elements from *ali* to their 5 prime end.

    The resulting :class:`pybedtool.BedTool` is sorted.

    :param alis: An iterable of :class:`pysam.AlignedSegment`
    :rtype: :class:`pybedtool.BedTool`

    :class:`pybedtool.BedTool` is an iterable:
    It can be iterated over in a for loop, or produce an iterator
    using :func:`iter` on it
    (see https://stackoverflow.com/a/28353158/1878788).
    """
    return BedTool(
        "\n".join(map(collapse_ali, alis)),
        from_string=True).sort(stream=True)

def main():
    """Main function of the program."""
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        "-b", "--bamfile",
        required=True,
        help="Sorted and indexed bam file containing the aligned reads.")
    parser.add_argument(
        "-o", "--out_bed",
        required=True,
        help="Bed file in which the results should be written.")
    args = parser.parse_args()
    with AlignmentFile(args.bamfile, "rb") as bamfile:
        collapse_and_sort(bamfile.fetch()).saveas(args.out_bed)
    return 0

if __name__ == "__main__":
    exit(main())

bam25prime/libbam25prime.pyx:

#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
# cython: language_level=3
"""
This library contains a cythonized version of a function to collapse aligned
segments to their 5' end.
"""

from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

def main():
    """Example use case."""
    print("No examples to show here.")
    return 0

if __name__ == "__main__":
    import sys
    sys.exit(main())

(Any advice or suggestions welcome: I'd like to improve my python package-making skills.)

For those interested, to use cython with pybedtools I had to define the following Extension (where glob is glob.glob and /home/bli/src/bedtools2-2.27.1/ is the location of the source code of bedtools).

    Extension(
        "bam25prime.libcollapsebed", ["bam25prime/libcollapsebed.pyx"],
        include_dirs=glob("/home/bli/src/bedtools2-2.27.1/src/utils/*"),
        language="c++"),


How to get the uniparc ids for all members of uniref cluster in batch request?

I found a way to batch the mapping from uniref cluster ids to the ids of all members. The url
https://www.uniprot.org/uniref/?query=UniRef100_A0A003&format=tab&columns=id,members
can be used to get the link from one uniref Cluster to it's member ids. Multiple uniref cluster can be fetched at once by connecting the ids with a +OR+ like so https://www.uniprot.org/uniref/?query=UniRef100_P99999+OR+UniRef100_P31946&format=tab&columns=id,members
Because the uniref cluster ids are not stable I ended up using the following syntax, where instead of querying the UniRef cluster by it's id, I query by one of it's members, namely the one that was used to create the current uniref cluster id: https://www.uniprot.org/uniref/?query=*&format=tab&limit=10&columns=id,members&fil=identity:1.0+AND+member:(P99999+OR+P31946)
In the notes section of the 'Retrieve / ID mapping' documentation of the uniprot they talk about fetching chunks up to 20.000 identifiers at a time. I am not sure if I use the service they are talking about though, so I'm keeping my requests at 100 identifiers per request and that seems to work.
One thing I noticed was that fetching results via the uniref endpoint is faster for small sets of identifiers (<20) since the fetch time of the batch request is similar for smaller and larger (<100) sets. So batching only makes sense for more than 20 identifiers at a time. (Disclaimer: I have only done very basic tests on this)
Once I have all the members of a cluster I can use the proteins/uniparc endpoint of the proteins api to fetch the uniparc data. For uniparc accessions (starting with the string 'UPI') I use the GET /uniparc/upi/{upi} endpoint. For the other accession, which should all be Uniprot accessions, I use the GET /uniparc/accession/{accession} endpoint.
This doesn't answer the question of a direct mapping between the uniref ids and uniparc member ids, but it is a workable solution.


nanopore - where to retrieve information from the basecaller used

I basecall separately with guppy. In the output folder specified by --save_path or -s there are a whole bunch of .log files.

$ ls -l *.log | head
-rw-r--r-- 1 tom tom 5242714 Dec  3 11:04 guppy_basecaller_log-2019-12-02_22-02-36.log
-rw-r--r-- 1 tom tom 5242718 Dec  3 11:06 guppy_basecaller_log-2019-12-02_22-04-38.log
-rw-r--r-- 1 tom tom 5242730 Dec  3 11:08 guppy_basecaller_log-2019-12-02_22-06-41.log
-rw-r--r-- 1 tom tom 5242879 Dec  3 11:11 guppy_basecaller_log-2019-12-02_22-08-50.log
-rw-r--r-- 1 tom tom 5242811 Dec  3 11:13 guppy_basecaller_log-2019-12-02_22-11-06.log
-rw-r--r-- 1 tom tom 5242817 Dec  3 11:15 guppy_basecaller_log-2019-12-02_22-13-28.log
-rw-r--r-- 1 tom tom 5242880 Dec  3 11:18 guppy_basecaller_log-2019-12-02_22-15-47.log
-rw-r--r-- 1 tom tom 5242816 Dec  3 11:20 guppy_basecaller_log-2019-12-02_22-18-09.log
-rw-r--r-- 1 tom tom 5242777 Dec  3 11:22 guppy_basecaller_log-2019-12-02_22-20-30.log
-rw-r--r-- 1 tom tom 5242742 Dec  3 11:25 guppy_basecaller_log-2019-12-02_22-22-47.log

The information you need for medaka is at the top of the first log file (alphabetically, and by creation date).

$ head -n 13 guppy_basecaller_log-2019-12-02_22-02-36.log
2019-12-02 22:02:36.255793 [guppy/message] ONT Guppy basecalling software version 3.4.1+ad4f8b9
config file:        /opt/ont/guppy/data/dna_r9.4.1_450bps_hac.cfg
model file:         /opt/ont/guppy/data/template_r9.4.1_450bps_hac.jsn
input path:         raw/minion
save path:          test_341/minion
chunk size:         1000
chunks per runner:  512
minimum qscore:     7
records per file:   4000
num basecallers:    1
gpu device:         cuda:0
kernel path:
runners per device: 4

If you ran live basecalling and you don't have those log files you could also try a recursive grep from the output folder, e.g.

grep -nr "^model file" .


What does presence/absence of black lines mean in IGV?

Black line represents an indel.  Read was lengthened to fit consensus.
No black line (blank), no read/sequence.  Read likely ended or was truncated.


What is the minimum number of individuals per group required for ADMIXTURE analysis?

The following paper provides an in-depth exploration of factors affecting power of admixture analyses under a number of different models. Although the authors state they "provide the open-source R software AdmixPower, to perform power and sample size analysis", in fact it seems the software is not available.

AdmixPower: Statistical Power and Sample Size Estimation for Mapping Genetic Loci in Admixed Populations

Tool availability:

https://research.cchmc.org/mershalab/Tools.html

Site asks for a password, apparently there is no option to sign up.


Bash scripting FastQC for multiple fastq files in multiple directories

multiqc kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the multiqc report and the actual fastqc html report to get a better idea of any error modes.

Going off of @Kubator's answer, I noticed that there was no command to run fastqc.

Here's a simple one-liner to run fastqc in parallel on all of your fastq files. The -j 25 uses 25 threads. Change 25 to however many threads you want/have for max speed.

#Run fastqc on everything in parallel.
> find ../reads/ -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' | parallel -j 25 --verbose
# copies all the fastqc files to directory ./
> find ../reads/ -name '*fastqc.*' | xargs -I '{}' mv '{}' ./

These files might be output by multiqc, anyway!

edit - now the command is safe for whitespace in file names. Please, no newlines in filenames though!


How can I find the chromosomal location of a list of genes?

Ok so I found two websites that could help. The one I found first is more intuitive to go through, but it looks like it could be outdated as the second website I found seems more recent so it probably has a more complete genome to use as the reference. (Sorry, I haven't worked personally with sunflowers and don't know for sure where you would go for the "official" reference genome, or if there even is an "official" one, but hopefully one of these websites will help.)

1) This is a website I found that looks like it is a website for the sunflower genome and includes gene annotation information. (The recent publications list only goes up till 2015 though so it may not be the most recent data)

https://www.sunflowergenome.org/annotations-data/

If you scroll down to the bottom of this page there is a text file detailing the descriptions of how the files are formatted.

For example, the linked to file gives the following:

  Ha412v1r1_CDS_v1.0.fasta.gz

  Description: Nucleotide CDS for each gene
  Header format: >TranscriptID GeneID Chromosome_start_end_strand
  Example: >Ha1_00044280-RA Ha1_00044280 Ha1_94590914-94591860_-

If you download this fasta file you could then create a python script that will read in the file and extract the headers. If each header is a string variable called "mytext", you could then extract the TranscriptID, GeneID, and chromosomal info and using mytext.split(" "). You could then split on "_" to further parse the chromosomal location info.  You could save the extracted information in a dictionary using the gene name as a key and parse through the original file you mentioned containing the list of genes, in order to get the chromosomal location information you desire.  The one caveat is that this is assuming that your gene names are the same as the geneIDs the above website uses.

I did not exhaustively examine the list of files and perhaps there is another one that would be a better one to parse through for your purpose, but the linked to text file gives more details on the files and their formats and presumably they could be parsed in a similar manner.

2) Double checking my first answer I found this article:

  Badouin H, Gouzy J, Grassa CJ, Murat F, Staton SE, Cottret L,
  Lelandais-Brière C, Owens GL, Carrère S, Mayjonade B, Legrand L. The
  sunflower genome provides insights into oil metabolism, flowering and
  Asterid evolution. Nature. 2017 Jun;546(7656):148.

Which links to a website as follows:

  "The integrative web interface Heliagene provides visualization, querying tools for data mining and network exploration for the community (https://www.heliagene.org)."

So it looks like this website also has reference genomes for sunflowers, although a lot of the annotation information seems to be locked unless you have a login. Like I mentioned previously, this website is a little less intuitive to go through but you might be able to find something helpful here if you think the other website isn't official or recent enough.


Two of those command line parameters should be accessible locations (see here). Could you please add the output of the following commands:
ls -l /data/fi1d18/Downloads/FASTQ_scATAC/ | head
and
ls -l refdata-cellranger-arc-GRCh38-2020-A-2.0.0 | head
I suspect that the refdata link should be an absolute path, rather than a relative path, i.e. something like this:
cellranger-atac count --id test --reference /opt/refdata-cellranger-arc-GRCh38-2020-A-2.0.0 --fastqs /data/fi1d18/Downloads/FASTQ_scATAC/ --chemistry ARC-v1


The sure why any software should be confined to only broad peaks, but there is two packages that should do what I thing you have in mind. Just to clarify, for smallRNAs peak detection is sometimes also referred to as loci or cluster detection:

ShortStack which is developed for alignment, annotation, and quantification of small RNAs. The de-novo cluster finding option as copied from the manual might do the trick:

DE-NOVO CLUSTER FINDING
    Unless options --locus or --locifile are used (see below), ShortStack
    will de-novo identify clusters of small RNA accumulation genome-wide.
    Cluster definition is simple: First, all regions containing at least one
    primary alignment are found where the maximum distance between the ends
    of the alignments is <= option --pad (default: 75). Second, if the
    number of alignments in the cluster is >= option --mincov (default:
    0.5rpm), the cluster is kept. The mincov threshold can also be specified
    in terms of reads per million by using a value such as 3.2rpm (which
    specifies the threshold to be 3.2 reads per million). Using a rpm
    threshold allows the sensitivity of cluster discovery to be normalized
    between libraries of different sizes. Alternatively, reads per million
    mapped (rpmm) can be specified: A --mincov of 1.2rpmm indicates 1.2
    reads per million mapped is the threshold. rpm is a fraction of total
    library size, while rpmm is a fraction of only the aligned & placed
    fraction of the library.

segmentSeq
is an R-package useing a Bayesian approach for small RNA locus detection. It's however really slow and not developed any further by the main author.


Detecting structural variants with MinION data

There was a Structural Variant breakout session at the London Calling conference this year. Unfortunately I didn't attend that session, but MinION community members have access to Constance Donnell's summary of that:

https://community.nanoporetech.com/posts/breakout-structural-varia

Here are my attempts at grabbing non-creative chunks from those notes:

Professor Wigard Kloosterman has been developing a bioinformatics pipeline called NanoSV for mapping genomic structural variants in patients with congenital abnormalities
Tomas Sesani at the University of Utah studied the precise copy number of duplicated genes in individual virus genes over time, and tracked the SNP within copy number variable regions of the genome [no public tool mentioned]
Dr. Sudha Rao from Genotypic Technology in India used the MinION sequencer with Sniffles for structural variation calling

ONT have said that the sessions/talks from London Calling 2017 will all be made publicly available at some point in the future.


Printing a bracketed list of numerical values is not typical for VCF—I'd guess it's not strictly allowed. A list should be printed only with the commas to separate values. Perhaps the solution is as simple as removing the square brackets.

cat input.vcf | perl -ne 's/CAF=\[([^\]]+)\]/CAF=$1/; print' > fixed.vcf

Replace cat with gunzip -c if your VCF file is gzip-compressed.

The Perl regular expression is

looking for the string CAF=[...]
grabbing the characters in ...
replacing the entire match with CAF=...

Here's a detailed breakdown of the Perl regular expression.

s/.../.../; — substitution syntax; replace the contents of the left block with the contents of the right block
[^\]]+ — one or more non-] characters; this is the trickiest bit; see this tutorial for more info on the [] syntax for character classes and the special ^ symbol for character class negation
(...) — the parentheses capture the matched pattern into the $1 variable
CAF=$1 — replaces everything matched in the first block with CAF= followed by the captured characters


How to apply upperquartile normalization on RSEM expected counts?

You can use the quantile function in base R to get the value of a particular quantile (e.g. 0.75 for the upper quartile). This can then be used as a factor for normalisation: divide the observed expression by this number.

## set up simulated data
data.mat <- floor(10^data.frame(sapply(1:5,function(x){rnorm(100)})));
colnames(data.mat) <- paste0("S",1:5);
data.mat[data.mat > 0] <- floor(data.mat[data.mat > 0]);

It's generally a good idea to first remove genes that have no counts in any samples. Rows can be removed using a logical selection based on the total expression in any given row.

## Only retain rows that have a positive total expression per gene
data.mat <- data.mat[rowSums(data.mat) > 0,];

Two methods of upper quartile determination are presented here, one which is based on the total counts, and one based on the expressed counts (i.e. excluding values with no expression). These both use the quantile function of R. The apply function calculates statistics that are based on a single specified dimension (or multiple dimensions) of the data. The first dimension of a 2D matrix in R corresponds to the row, while the second corresponds to the column. By using 2 in the apply function, data are collected along the second dimension (i.e. the column), and statistics generated from that:

data.quantileAll <- apply(data.mat, 2, function(x){quantile(x, 0.75)});
data.quantileExpressed <-
  apply(data.mat, 2, function(x){quantile(x[x>0], 0.75)});

Normalisation is applied by dividing the values by the column-based upper quartile statistic. R will repeat values along the last dimension first (which in this case is the column), so the matrix needs to be transposed first (with columns as genes) to get the arithmetic working properly, with a final transposition to return the orientation back to the original direction (with columns as samples).

data.norm <- t(t(data.mat) / data.quantileExpressed);


Single cell RNAseq cell cluster (true cluster or sub cluster)

I wouldn't be surprised if the 15 cell types that you are expecting were mostly characterized via "protein-based" methods; antibody staining, flow cytometry using fluorophore coupled antibodies, ... The data that scRNA-seq is RNA based and should not be expected to provide the same information as protein based assays. Moreover, the scRNA-seq data, especially those from microfluidic systems like 10x, are sparse and might not be enough to resolve closely related cell types/states such as those of T cells.
On top of the aforementioned, evaluating clusters or clustering stability is a difficult task and in my experience, quite some cluster evaluation metrics cannot be easily applied to single cell data simply because of the shear size of the data. Having asked your question and looked for answers (and I am still doing that), I am doing the following:
i) Try to use transcriptomic markers as much as possible. Even then a lack of a marker would not mean much in terms of scRNA-seq, it might very well not be detected simply due to sensitivity
ii) Use clustree do select the most plausible number of clusters (resolution parameter in Seurat). The package is compatible with Seurat and some other scRNA-seq packages.
iii) Use silhouette width as a metric for clustering. This is computationally expensive, however, should be fine for your 5000 cells. For larger cell numbers, I use another package that approximates this (will add a link but first I have go through my scripts for the name).
iv) Check if batch effects are accounting for what you call "sub-clusters". For example clusters 0, 14, 24 and 25 above might correspond to different samples processed in different days.


Getting a "system is computationally singular" error in sleuth

This happen when the variables (strain +batch) create a design matrix like this:

batch strain
1 1 #
1 1 #
1 2
2 2
3 3
4 3
...
16 72

Which means that some of the covariates are not linearly independent (ie batch 1 and strain 1), all the strain 1 is in batch 1.

You can correct for batch effects, but not in this design of the linear model (if you want to take into account the strain). You could do one batch more with those 15 strains that are in a single batch (if they are in different batch between them) that way you would get an independent design.

Many strains in a single batch are in the same batch. You need to increase the number of samples (recommended anyway due to the low number of samples per strain) to avoid this problem.

There are a lot of related question in Bioconductor support forum, from where I expanded an answer.


How can I compare structures of proteins responsible for dehydration tolerance in plants?

Have you explored in PDB?. If you already have a starting sequence, you can use it as input to search with phmmer using PDB as database. Once you have a curated set to analyse, you can visualize and align the structures with VMD, for example.


Get Gene Expression Matrix from GEOquery


I wrote a command-line k-mer counter called kmer-counter that will output results in a form that your Python script can consume: https://github.com/alexpreynolds/kmer-counter

You can grab, build and install it like so:

$ git clone https://github.com/alexpreynolds/kmer-counter.git
$ cd kmer-counter
$ make
$ cp kmer-counter /usr/local/bin

Once the binary is in your path, you might use it in Python like so:

k = 6
fastaFile = '/path/to/some/seqs.fa'
kmerCmd = 'kmer-counter --fasta --k=%d %s' % (k, fastaFile)
try:
    output = subprocess.check_output(kmerCmd, shell=True)
    result = {}
    for line in output.splitlines():
        (header, counts) = line.strip().split('\t')
        header = header[1:]
        kmers = dict((k,int(v)) for (k,v) in [d.split(':') for d in counts.split(' ')])
        result[header] = kmers
    sys.stdout.write("%s" % (str(result)))
except subprocess.CalledProcessError as error:
    sys.stderr.write("%s" % (str(error)))

Given example FASTA like this:

>foo
TTAACG
>bar
GTGGAAGTTCTTAGGGCATGGCAAAGAGTCAGAATTTGAC

For k=6, you would get an iterable Python dictionary like this:

{'foo': {'TTAACG': 1, 'CGTTAA': 1}, 'bar': {'GTTCTT': 1, 'AGAACT': 1, 'GAGTCA': 1, 'ATGGCA': 1, 'GAACTT': 1, 'ATTCTG': 1, 'CTAAGA': 1, 'CTTCCA': 1, 'ATTTGA': 1, 'GGAAGT': 1, 'AGGGCA': 1, 'CCTAAG': 1, 'CTCTTT': 1, 'AATTTG': 1, 'TCTGAC': 1, 'TTTGCC': 1, 'CTTAGG': 1, 'TTTGAC': 1, 'GAAGTT': 1, 'CCCTAA': 1, 'AGAATT': 1, 'AGTCAG': 1, 'CTGACT': 1, 'TCTTAG': 1, 'CGTTAA': 1, 'GTGGAA': 1, 'TGCCAT': 1, 'ACTCTT': 1, 'GGGCAT': 1, 'TTAGGG': 1, 'CTTTGC': 1, 'TGGAAG': 1, 'GACTCT': 1, 'CATGCC': 1, 'GCAAAG': 1, 'AAATTC': 1, 'GTCAAA': 1, 'TGACTC': 1, 'TAGGGC': 1, 'AAGTTC': 1, 'ATGCCC': 1, 'TCAAAT': 1, 'CAAAGA': 1, 'AACTTC': 1, 'GTCAGA': 1, 'CAAATT': 1, 'TAAGAA': 1, 'CATGGC': 1, 'AAGAAC': 1, 'AAGAGT': 1, 'TCTTTG': 1, 'TTCCAC': 1, 'TGGCAA': 1, 'GGCAAA': 1, 'AGTTCT': 1, 'AGAGTC': 1, 'TCAGAA': 1, 'GAATTT': 1, 'AAAGAG': 1, 'TGCCCT': 1, 'CCATGC': 1, 'GGCATG': 1, 'TTGCCA': 1, 'CAGAAT': 1, 'AATTCT': 1, 'GCATGG': 1, 'ACTTCC': 1, 'TTCTTA': 1, 'GCCATG': 1, 'GCCCTA': 1, 'TTCTGA': 1}}

You can use standard Python calls to manipulate this dictionary object and get sums of counts per record, for sequence, etc. which seems to answer your question. Please feel free to clarify what you're looking for if this object representation is not clear.

For a fully-fleshed out demonstration, see: https://github.com/alexpreynolds/kmer-counter/blob/master/test/kmer-test.py


Where can I find gene expression data on one of the cell lines in NCI-60?

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

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


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


Parallel processing of scripts that use obitools

Your only option here is to first split the fastq files into small chunks, process each and then merge them back together. Alternatively, find something else that runs faster.


Concerning the current outbreak, using bioinformatics

My personal approach to something like this is to first try to reproduce existing work in order to learn to use relevant tools and understand the concepts involved. Then as you get more comfortable with the subject matter, you may spot an opportunity to add value.

If you'd like to get started with the phylogenomics side of things, you can use this public workspace developed by the viral genomics group at the Broad Institute, which includes viral genome sequences from NCBI SRA as well as some workflows for doing genome assembly and quality control: https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19


How many models should be generated by rosetta for denovo protein structure prediction?

Unfortunately you do in general increase your chance of success by going to the 50k range of models or higher. Projects like Rosetta@home exist because even large research groups can't access as much compute power as they'd like. In this case the long time arises from the combinatorial explosion of possible fragment combinations that need to be explored.

Success is rather protein dependent though - sometimes the right structure falls out after only a few models, other times your fragment library does not cover the protein adequately and you could run Rosetta forever without getting a good structure. Have a look at your 1,500 models and see how similar they look.

382 residues is rather long for Rosetta. Another factor to consider is whether you are using predicted contacts as part of the folding process (e.g. in the PconsFold protocol), which will greatly improve the chances of getting a good model and reduce the number of models you need. In fact the PconsFold group switched to CONFOLD for PconsFold2, I think mainly because it was much faster, so you are not alone in finding Rosetta runs long.


counting individuals with the present gene and absent gene for each population

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

Here an example for the KS2W files:

myFiles <- dir(pattern="KS2W")

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

fileNames <- list()

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

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

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

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

output <- cbind(present, zeros)

rownames(output) <- merged[,1]

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

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

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

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

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


Understanding some of the computational bottlenecks of Covid-19 research

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.


Pymol: color spheres with different colors

Spheres
So how PyMOL stores how it shows atoms is that each atom has a colour and a represention, which is a binary number telling what to show. Therefore the colour for spheres is the same as the colour for sticks. The only way round this is to copy the object (create copycat, original and you set copycat differently)
To set the atom you do
select iliketheschaps, resi 200‐221 and name CA+CB
color tv_green, iliketheschaps
show spheres, iliketheschaps

The object iliketheschaps is a selection. Do make sure the range is a hyphen-minus (character next to 0) as opposed to an en-dash (–) which Word subs and on a mac you get with cmd+-
non-PyMOL colours
Now you asked for "gold". This is not a colour in the menu —famous for the colour limon. Gold is a CSS colour name, but it won't work.
You will have to colour by hex code. Say you want "gold" you need to google the colour name and hex. #FFD700.
color 0xFFD700, iliketheschaps

The Fashion Police have been alerted

This paragraph is in jest.

However, as you may know bioinformaticians are very fashion savvy for colours. Mustard was the in-colour from last year, so you really ought to use mustard (0xffdb58). This year is all about nectarine (0xffa168) (which is not salmon —this ain't 2015!). Coral (0xff7f50) and teal (0x008080) were from previous years and are still are going strong!


Why are inversions defined as the reverse complement and not just the reverse of the reference?

If an inversion were just reversed then there would be 3' -> 3' bonds and 5' -> 5' bonds. That's why inversions are reverse complemented, you then maintain the normal 5'->3' direction. This is nicely covered over on our biology sister site, which also includes this image:


How to calculate overall reference coverage with MUMmer?

I believe all you need to do is to run dnadiff from MUMmer. That will run a comparison and output a number of useful metrics.


Parsing specific features from Genbank by label?

The /label=5'ITR is called a qualifier. You can look for those qualifiers:

for feature in reference.features:
    for k, v in feature.qualifiers.items():
        if k == "label" and v == ["5' ITR"]:  # or "5' ITR" in v
            print(feature.location)


How to get UCSC Genome Browser to show the actual sequence?

If you zoom in at one specific region, it shows the nucleotide chain on the third line (below the chromosome number). It can't be there by default as the size of default region is more than 100,000bp. You can use the scale on the top of the page to see the size of the region.


plink: --update-name vs. editing the BIM

I would reccomend using bcftools on the original vcf files before you convert them to plink, to fill in missing IDs using the command:
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf

This means you won't have any duplicate ID names and plink should be happy about that. It will preserve the positional information in the .bim file as well.


how to count the number of each allele of CNVs?

Here's a pure python solution to process your example. It should be straightforward to loop over all the rows of your vcf files

my_line='16 32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3 '
my_line=my_line.strip() # remove trailing whitespace
my_fields=my_line.split()
names=my_fields[4].split(',') # extract CNV names
names=[x[1:-1] for x in names] # remove '<>' symbols
print('\t'.join(names))
for sample in my_fields[9:]:
    genotypes=[0]*(len(names)+1) # initialize array of zero counts, including reference
    for allele in sample.split('|'):
        genotypes[int(allele)]+=1
    print('\t'.join([str(x) for x in genotypes[1:]])) # print all columns except reference


Where are .motif files from homer knownResults?

I have provided a solution in Biostars:

If you edit the file findKnownMotifs.pl and move the call to printMotif to come before testing $svgFlag you will find that .motif files are printed again.

In the following I have commented out the call from where it need not be, and introduced it as the first line, right before the if block:

   printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
    if ($svgFlag) {
        my $str = HomerSVGLogo::motif2svg($_,$bitFlag,$rnaFlag2);
            open SVG, ">$directory/$knownDir/$logoFile.svg";
            print SVG $str;
            close SVG;
            print MAIN "$str\n";
    } else {
  #printMotif($_, $directory . "/$knownDir/$matrixFile"); # this should be done regardless of $svgFlag
            `profile2seq.pl "$directory/$knownDir/$matrixFile" 100 > "$tmpFile"`;
            my $width = $widthFactor * $_->{'len'};
            `seqlogo -a $bits -f "$tmpFile" -F PNG -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width`;
  `seqlogo -a $bits -f "$tmpFile" -F PDF -S -c -o "$directory/$knownDir/$logoFile" -h $height -w $width 2> /dev/null`;
  print MAIN "<IMG src=\"$knownDir/$logoFile.png\"/>";
}

I expect this change was mistakenly introduced while refactoring the code to support svg as alternate output format for logos.

I have reported this to the author of this code.


How to combine two Genome-wide Association Study (GWAS)?

The first way is to combine two samples and re-run GWAS. A second way is to do a meta-analysis using GWAS results.


How can I extract the longest N isoforms per gene from a fasta file?

I know you are asking for a python solution. But maybe you are interested in other ways.

You can use samtools faidx to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.

1. Index fasta

samtools faidx input.fa

2. Extract IDs of the n longest transcripts per gene

sed 's/_.\+//g' input.fa.fai | sort -u |
    parallel -k 'grep {} -m 3 <(sort -k2,2 -r --sort g input.fa.fai) |
        cut -f1' > longestID.txt

The sed command replaces everything from _ onwards in each line, leaving behind only the gene name. sort -u is used to get a list of distinct gene names. For each of theses gene names we grep for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with cut. Adopt the -m parameter if you want more or fewer transcripts.

3. Get sequences

samtools faidx input.fa -r longestID.txt > output.fa


Interpreting Intergrative Genomic Viewer (IGV)

You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad institute software engineers replying here.

Your next question is likely why the color of the reads and the orientation of the genes are opposite of each other. This is because in modern directional protocols, read #2 in a pair matches the orientation of the gene from which it arose.


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

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

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

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

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

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

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

Hope that it answers your question.


Why is cell count stated relative to animal size?

The typical goal of scRNA-seq is to assess the possible cell-types and states that exist in an organism/tissue/etc. Given that, it can be useful to know what portion of the possible cells have been assayed, since that sets bounds on the number and size of cell populations not assayed. If the goal of a study is to assess the variability in a single cell-type, then you're correct that mentioning "0.8% of the population was sequenced" isn't actually informative. Note that "in a single cell-type" is important there, since if one is assaying multiple cell-types then it's important that enough of each are assayed given their prevalence in the relevant population.

Practically speaking, numbers like this are only ever written in papers if they look impressive, so no one is going to write about sequencing 0.000001% of the cells in an organism, since the reader/reviewer will then start to question the relevance of the results.


DMPfold error while running protein structure prediction

The comment by Devon Ryan is correct - looks like your blastpgp command is not working properly. For a protein of that length it's more likely to be a problem with the BLAST executable running on your system than a memory issue.

Try running the blastpgp command (https://github.com/psipred/DMPfold/blob/master/bin/runpsipredandsolvwithdb#L26) on its own and see what happens:

~/path/to/blastpgp -a 12 -b 0 -j 3 -h 0.001 -d yourfile.a3m -i yourfile.fasta


How do you currently track bugs in python using freely available software/websites?

Use IPython and activate the trace in the Python debugger as follows
%pdb 1

or
%pdb on

to turn this off %pdb 0 or %pdf off
I do recommend using a formal IDE such as PyCharm or Visual Studio Code. The professional version of PyCharm is highly rated and within that use breakpoints . In Python3.7 there is supposed to be breakpoint() so you may not need an IDE, i.e. you set the breakpoint in your code using breakpoint(). In older versions of Python you can import this function via, import pdb; pdb.set_trace()
Finally check out PDB++ which is found here
Its important to note that experienced programmers often debug by simply inserting print() throughout their code and don't bother with an IDE


How to predict protein binding to ligands?

If you are training a machine learning algorithm you probably want to train against things you know to be true, rather than making predictions from the sequence and training against the predictions.

In this case you will want to use protein structures with the relevant ligands bound as your data. If you go to the RCSB PDB advanced search you can search for entries with any ligand present, specific ligands present, protein and DNA present, etc. Then you can extract the FASTA sequences for these entries.

See also databases like MetalPDB for specific ligand types, and one of our papers for an example of creating a dataset of metal binding proteins.

Features to use would be the raw sequence as a one-hot encoding on amino acid types. You could add other sequence level features such as secondary structure predictions, disorder predictions or physico-chemical properties. However a decent network would be able to predict much of this if it needs to. You could create alignments and put in a sequence profile or covariation data, which would be giving the network more information.


CRISPR Sequence Finder and Database Download

CRT (CRISPR Recognition Tool) http://www.room220.com/crt/

Crisprs Finder Online Tool https://crispr.i2bc.paris-saclay.fr/Server/

Piler-CR http://www.drive5.com/pilercr/

(Omic Tools List) https://omictools.com/crispr-detection-category


How do you download stuff from NCBI fast?

Answers to this question can be found on biostars. To give you a short summary:

Faster way might be to use the parallel-fastq-dump, as suggested in this answer. I never tested that tool though, my own experience is that prefetch is more stable than fastq-dump command.
In that same post it seems that the Japanese mirror seems faster, don't know if that is still valid.
The reason this is never optimized might be that downloading from SRA is not a task you'll have to do repeatedly, as described here in this post.


Downloading URL using AWK for fixed fields

awk '{print "wget " $9 " -O " $1 "_" $2 ".fastq.gz"}' tmp

will give you

wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz -O b_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz -O c_k.fastq.gz

If you add a | bash to the command it will then execute all the wget commands.


How can I use Nanopore reads to close gaps or resolve repeats in a short-read assembly?

You might want to look into Unicycler (manuscript with more information can be found here); even though it is supposed to be used with bacterial genomes only, it might work well with a small genome such as a mitochondrion's.

Beware that if you happen to have very long reads, you might end up with an assembly with multiple copies of the circular genome: you might want to look into circlator then.


Grouping a set of genes in a ggplot object

If you had the classification into oncogene/TSG for each gene in a column (let's call it gene_class) in lmuts, you could split both groups into facets within the plot by adding:

+ facet_wrap(~gene_class, ncol=1)


How to remove sequences from a fasta file using a sequence ID list which contains a space within the id?

Install (pip) biopython if u don't have it.
from Bio import SeqIO
with open('/path/ids.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('/path/output.fa', 'w') as fout:
    for record in SeqIO.parse('/path/sequences.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')


c++ program to find length string an char*

I think you're confusing C++ and C strings, which are very different. When you have a char *, that's a C-style "string". It's really just a bunch of bytes with a zero byte at the end, demarcating the end of the string. To find the length of a C string you include <cstring> (or <string.h> in C) and call strlen, like this:

/** Print every character in the string */
void print_chars_in_string(char* word) {
    for (size_t i = 0; i < strlen(word); ++i) {
        putchar(word[i]);
    }
}

strlen is actually pretty simple though, you can write your own version as an academic exercise like this:

size_t deluxe_bespoke_strlen(char* str) {
    size_t length = 0;

    while (*str++)
        ++length;

    return length;
}

Since the while-loop condition will be true until you hit the null terminator (the zero byte at the end), this is how you would implement the strlen function without a for-loop, otherwise you would need the strlen function for the strlen function.

In C++, strings are actual objects, with methods you can call. To get the length of a C++ string simply do this:

std::string my_name = "Jose";
std::cout << my_name.length() << std::endl;

If you have a C string, you can convert it to a C++ string with the built-in constructor:

const char* my_name = "Jose";
std::string my_name_cppstring = std::string(my_name);

You can also get the C version of a C++ string by calling the c_str() member function like this:

std::string sequence = "GATTACA";
size_t length = strlen(sequence.c_str());

My recommendation is to not mix C and C++. C is my favorite programming language, but that doesn't mean everyone should write everything in C. Since you're using C++ I'm guessing that's your preference, so I would suggest using C++ strings instead of C strings, just remember to include <string>.

As for your specific question, the compiler is giving you the correct message: str is not a member of std. You don't add the namespace to variables. This is how that works:

#include <cstdlib>
#include <iostream>

namespace CustomDeluxeBespokeStandardNamespace
{
    size_t custom_bespoke_strlen(char* str)
    {
        size_t length = 0;

        while (*str++)
            ++length;

        return length;
    }
}

int main(int argc, char *argv[])
{
    if (argc == 1) {
        std::cerr << "No input(s)." << std::endl;
        std::exit(EXIT_FAILURE);
    }

    while (*++argv) {
        std::cout << *argv << ": " << CustomDeluxeBespokeStandardNamespace::custom_bespoke_strlen(*argv) << std::endl;
    }
}

So if you run it like this:

testprogram horse house mouse cat

You should get this output:

horse: 5
house: 5
mouse: 5
cat: 3

So to recap, you were adding the namespace to your variables, and the compiler was giving you an error because it correctly could not find anything by that name in the std namespace. Even if it had though, C style strings do not have member functions. There are no such things as objects in C, at least not in the sense that you can define member functions for an object, so that was the second error.

As an aside, you will note all my examples using strlen both custom and actual return size_t. This is not an accident; in both C and C++ the return type of a call to strlen is size_t. You have to be careful because in your code you compare an int to a size_t. Comparing a signed value to an unsigned value is Not Good because it will work only if the signed value can be faithfully represented after a conversion to an unsigned value. If it can't though, you're in for a bad time.


Are phylogenetic tree construction algorithms any different than general clustering algorithms?

The goal of a phylogeny is to estimate the "expected" number of mutations between all taxa in the analysis and their hypothetical common ancestors. A cluster-analysis will only identify the "observed" mutations and "expected" and "observed" mutations can be majorly different due to the major artefact of reversion mutation. This is particularly true of nucleotide phylogenies.

The key difference between clustering algorithms based on a "uncorrected" distance matrix and phylogeny is the latter is based on an explicit model to accommodate reversion mutations. The real problem is that there are only 4 bases, so randomly there is 1/4 chance that a mutation at a given position will mutate back to the original, e.g. A->C->A. Observed mutational differences = 0, expected (real) mutational differences = 2. What phylogeny is concerned about is reconstructing that "2". The problem is significant because almost every gene has regions of rapid mutations and regions of low mutations.

The primary way of doing this is via the Jukes-Cantor correction and is foundational in all nucleotide trees expect "p-distances". If nucleotide divergence is lower than 75% then it is possible to estimate the expected number of mutations using the observed via the JC correction. Moreover, when combined with a method of estimating rate variation within a gene, (usually the discrete gamma distribution), JC correction is very effective at recovering the "true tree". This is because fast evolving homologous sites are grouped together - big correction via JC, slow evolving sites are grouped together - small correction via JC. Other approaches to improving the JC correction are via identifying the bias between purine to purine mutations and pyrimidine to pyrimidine and purine to pyrimidine mutations.

The importance of JC was demonstrated by simulation studies of 4 taxa ((a,b), (c,d)), if two strains evolved very quickly when there sister lineages evolve slowly a clustering algorithm will report the fast strains are a sister group, i.e. ((b, c),(a, d)). If the JC method is implemented via maximum likelihood (or Bayesian) then it correctly recovers the true tree ((a,b), (c,d)). The artefact is known as long branch attraction.

Clustering algorithms based on a distance matrix that implements the JC correction, tend to perform poorly for long branch-attraction artefacts. This doesn't mean the correction is useless, just not particularly powerful. The issue is that distance matrix methods don't "stick" to the model and the clustering will introduce a layer of inaccuracy. Normally presentation of a "corrected" distance matrix combined with a clustering algorithm will require bootstrapping (resampling with replacement) to assess whether a given cluster is supported. Parameterised distance matrices, using neighbour-joining clustering in in-conjunction with a bootstrap are considered okay.

The @user172818 mentioned parsimony and this method is considered less reliable because it can't implement a JC correction. IMO it is possible weighted parsimony could make a "come-back", but it would be really complicated to implement a biological weighting method and would require extensive, independent calculations.


Finding transposable elements using RepeatMasker

You need to model the repeats in your de novo genome.

See:
http://www.repeatmasker.org/RepeatModeler/ and
https://www.biostars.org/p/154290/ (from previous answer)

You first build a database named "name_of_your_database" (insert whatever you want to call it) and run it on the fasta file of your genome (e.g. your_genome.fasta).

Make you find the right paths for the software.

## build the database
BuildDatabase -engine ncbi -name name_of_your_database your_genome.fasta

## run the modeler (in this case using ncbi)
RepeatModeler -database name_of_your_database -engine ncbi


How do I add external labels and connecting lines to a Ternary Plot in R?

Hmmm I'm not familiar with AddtoTernary function, I can see that you've placed a function within the "AddtoTernary" function. The traditional ggplot2 way of doing this would simply be controlled by the theme

ggtern(data=x,aes(x1,x2,x3)) +
  geom_point(fill="red",shape=21,size=3) +
  theme_tern_bw()

This will give an arrow for each axis, but additional code is required for different colours (I'd have to look up my past work).

I think what you might want is a colour key? In this case ggtern will have a "key" function because ggplot2 has one.Here, your colouring is based on species not on a colour gradient as to how far to the corner of a graph you are.

Alternatively you can easily put names besides each point in ggplot2, but again this is traditional style code.

Also you can add different shapes to the graph, which might be useful instead of colours.

Arrows and names together peronsally I've not seen on a ggplot2 graphic and I'd have not idea how to do it.

FYI I don't think we need all the data just a few points fine, although it shows why you prefer arrows - personally I would use colour coding. In "traditional" ggplot2 coding the graphing is built up in layers, so if additional functionality is required it is simply appended to the exisiting command. Thus it is quite easy to make modifications.


How to retrieve InChI key for KEGG compound?

I'm not aware of a functionality in the KEGG service that would allow you to do that directly. I believe your solution (using e.g. ChEBI) is the correct one.

You may have other useful service such as unichem, chemspider (you would need a login though for that one) could help as well.

Disclaimer: I'm the bioservices main author


Get residue position from Pymol Selection

Ok, after some more searching I found the answer.

You have to use the cmd.iterate function to iterate over the default selection ("sele", first parameter). With the second parameter you can push the variables you need into an array for later access. In this case I wanted the residue number so I used the variable 'resv'. A full list of the variables exposed by the iterate function can be found on the pymolwiki

Note that you can use a function call as the second parameter and if you need a more complex selection than 'sele' you can use selector.process in the iterate function.

Here is my final working function. Note that I used a set for storing the residue ids, since the iterate function technically iterates over the atoms and so the same residue will likely appear multiple times in your selection.

def show_selection_position_data():
    residue_data = {
            "header": ['variation', 'rmsf']
            , 1: ['high', 0.1]
            , 2: ['high', 0.2]
            , 3: ['low', 0.1]
            , 4: ['high', 0.05]
            , 5: ['middle', 0.2]
            , 6: ['high', 0.1]
            , 7: ['high', 0.1]
        }
    selections = cmd.get_names('selections')

    if "sele" not in selections:
        print("No selection")
        return

    userSelection = "sele"

    stored.residues = set()
    cmd.iterate(selector.process(userSelection), 'stored.residues.add(resv)')

    for residue in stored.residues:
        if residue in residue_data:
            print(residue_data[residue])
        else:
            print("{}: no data".format(residue))

    return


Is the optional SAM NM field strictly computable from the MD and CIGAR?

Assuming both the MD and CIGAR are present and correct, then yes, you can parse both to get the edit distance (NM auxiliary tag). One big caveat to this is that there's a reason that the samtools calmd command exists, since it's historically been the case that not all aligners have output correct MD strings. It's rare for the CIGAR string to be wrong and that'd be more of a catastrophic error on the part of an aligner. For what it's worth, if the NM auxiliary is absent on a given alignment but present on others produced by the same aligner then it's fair to assume NM:i:0 for a given alignment by default (many aligners only produce NM:i:XXX if the edit distance is at least 1).


Why would someone use a CRAM instead of a BAM?

Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings.
If you need to guarantee that any random obscure downstream program will be able to handle it.

Uptake of CRAM has been pretty slow. Java programs using htsjdk (e.g., picard, IGV and GATK) have only relatively recently added support for CRAM. If you need to use an old version of those for some very odd reason then CRAM may not be supported.

There are a lot of programs written in python that use pysam to open BAM files and these should, theoretically, support CRAM. The issue is that some of the functions may fail and one can't assume that authors will have always written the code needed to handle this. I'll use deepTools as an example, since I'm one of its developers. One of the things about CRAM files is that they (by default) are made such that you require a reference genome in order to construct the sequence field in each alignment. This works fine if you're using a standard genome (htslib, via pysam, can fetch many standard genomes from the web automatically), but if you're not, then you need to specify a fasta file to use for decompression. Every tool, then, needs to add an option for this. With pysam 0.14 and htslib 1.7 this can be circumvented by not decompressing the sequence, but behavior has to be explicitly requested.

Another issue is that many tools will use features from the file index, such as the .mapped accessor, to get the number of mapped reads in a file. CRAM files contain very very little information, so this then fails. Consequently, tool authors need to check for CRAM files and both derive and propagate this information through their functions if it's needed. This can be a time-consuming task (e.g., it took me a couple days to get this implemented in deepTools). Relatedly, samtools idxstats is useless on CRAM files, since there are no statistics stored in the index.

That having been said, it's likely that CRAMs slowly gaining acceptance will eventually make it the standard. It's already a convenient archival format, it's just a matter of time before users can assume that most analysis programs are written to handle it.


tl;dr

Spliced aligners just make many small matches along the read pair or single read to determine splice junctions and output an alignment for the whole read, as long as it matches somewhere in the genome.

mRNA and translation

In the cell, after transcription, introns are spliced out of RNAs to become mature RNA. The protein-coding sequence does have sequence that is chimeric from cis-located portions of the same chromosome that do not actually exist.

              |<<exon 1>>|  intron    |<<exon 2>|
Genome: ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

mRNA  : ----ABCDEFGHIJKLYZ123456789aaaaaaaaaaa

As you said, if we only even consider small stretches of DNA, the sequence JKLYZ does not exist in the genome, only in the mRNA. So, we need a special way to map mRNA-seq reads to the genome that correctly splits at splice junctions.

How do we map mRNA-seq reads to a genome?

As a result, parts of the mRNA-seq reads will align well to one part of the genome, but only for part of the read. This requires a different alignment algorithm that considers multiple local high-score alignments for different parts of the read.

For example, see how a genomic aligner versus a splice aligner might treat the same mRNA-seq read on our example genome.

                         |<<exon 1>>|  intron    |<<exon 2>|
Genome:            ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

spliced aligner  : ------ABCDEFGHIJKL------------YZ123456789
genomic aligner  : ------ABCDEFGHIJKL

There is a possibility that the genomic aligner will not recognize that the second half of the read is an exon due to the large intron gap, and it may clip the alignment. Spliced aligners just take this into account and output several alignments for the same read regardless of gap size. The STAR algorithm, for example, uses many small matching pieces of a single read to determine where the splice junctions are. HiSat2 uses small approximate matches using a FM-index to find seeds for spliced alignment.


How to represent a deletion at position 1 in a VCF file?

From the latest VCF specs (page 8):

REF - reference base(s): Each base must be one of A,C,G,T,N (case insensitive). Multiple bases are permitted. The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event occurs at position 1 on the contig in which case it must include the base after the event;

So I think you're right


A quick look at your link tells me the SRR numbers run from SRR837819 to SRR837856. You can use fastq-dump from the sratoolkit, and make a for loop around it in bash.

Something like this should work:

for (( i = 19; i <= 56; i++ ))
  do
  fastq-dump --accession SRR8378$i
done

After reading Devon Ryan's answer, I realize that you asked for SRA files instead of fastq. This can also be done with prefetch of the sratoolkit:

for (( i = 19; i <= 56; i++ ))
  do
  prefetch SRR8378$i
done


efficient counting of dinucleotides/trinucleotides on fastq reads?

Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:

zcat file.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

I tested this on a 2.6G fastq file (exome data) and it took:

real    1m18.987s
user    1m41.392s
sys     0m10.684s

On a 66G WGS sequence file, it took:

real    34m6.771s
user    40m16.452s
sys     4m10.176s

The extra grep CGC isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the grep:

zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
    sort -rnk2

Following the OP's suggestion, I also tested with pigz and, it turns out, pigz can fly!

time unpigz -c exome.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

real    0m44.012s
user    1m19.544s
sys     0m14.900s

And for the 66G WGS:

real    17m29.935s
user    29m23.816s
sys     5m17.856s


What happens if a major bug is discovered in a bioinformatic package that has been used in published literature?

I prefer to treat software tools and computers in a similar fashion to laboratory equipment, and in some sense biology in general. Biologists are used to unexpected things happening in their experiments, and it's not uncommon for a new discovery to change the way that people look at something. Things break down, cells die off quicker on a Wednesday afternoon, results are inconsistent, and that third reviewer keeps on about doing that thing that's been done a hundred times before without anything surprising happening (just not this time). It's a good idea to record as much as can be thought of that might influence an experiment, and for software that includes any input data or command line options, and especially software version numbers.

In this sense, a discovered software bug can be treated as a new discovery of how the world works. If the discovery is made public, and other people consider that it's important enough, then some people might revisit old research to see if it changes things.

Of course, the nice thing about software is that bugs can be reported back to the creators of programs, and possibly fixed, resulting in an improved version of the software at a later date. If the bug itself doesn't spark interest and the program gets fixed anyway, people unknowingly use newer versions, and there might be a bit more confusion and discussion about why results don't match similar studies carried out before the software change.

If you want a bit of an idea of the biological equivalent of a major software bug, have a look at the barcode index switching issue, or the cell line contamination issue.


programmatic secondary structure prediction for >36-mer DNA oligonucleotides

If you really insisted on using Primer3 you could perform a 36bp sliding window analysis. You would submit your primers in 36bp 'chunks' with a step size of say 3 bp. You would need to "circularise" your primer within the sliding window to assess the "beginning" and the "end" and remove "artifactual" hairpins. Its doable via a script and you would automate the web upload and data extraction via a "digger/spider" or else just use a download equivalent (I've used a standalone). There would be a "blindspot" using this approach, which you could overcome if you wanted.

So summary, its complicated and I think your suggestion of using UNAFold is a better idea.


How to decide number of neighbors and resolution for Louvain clustering

A general rule of thumb is that in order to improve the variance $n$ times you need $n^2$ neighbours. This is only applicable if you consider the $n^2$ nearest neighbours of a cell to be biologically identical (i.e. "similar enough"); if your data includes 10 types of cells with 10 cells each, then using the 20 nearest neighbours for smoothing will obscure the data.

As far as I know, there is no single best answer to this question. I would suggest trying different numbers and sticking to what agrees more with the biology of the dataset.


Are there open-sourced graph-based variant callers?

Yes, there are. There are some suggestions in the comments (VG, WhatsHap, GraphAligner, Minigraph).
However, to be clear, the current default variant calling algorithm used in the GATK (the haplotype caller) is a graph-based algorithm. It constructs a De Bruijn graph from mapped reads for local assembly around variants. Other tools that also follow this strategy include Platypus and freebayes.
You might also consider multi-color genome assemblers. For example, Cortex and Corticall use colored De Bruijn graphs to aggregate reads from multiple genomes and call variants.


Albacore wheel not supported on this platform (Mac OS X El Capitan)

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

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

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

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

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

# activate the virtual environment
pipenv shell

read_fast5_basecaller.py -h

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