How GFF3 attributes (9th column) varies from one gene prediction algorithm to another GFF3 files are in tabular format with 9 fields per line, separated by tabs. The first 8 fields share almost same data structure, but the 9th field varies a lot depending on feature type and gene prediction algorithm.

Presently I am trying to build my own GFF module that will parse a GFF file and return one structured data object (a python dictionary structure). For that, I need a good understanding of the 9th field of the GFF format.

So can any one please help me by providing various GFF files?

The first place to start is the GFF3 specification. This is the official word on what is and is not allowed in a GFF3 file. For example, users can define arbitrary attribute keys, so long as they do not begin with an uppercase letter (these are reserved for "official" use).

But your question doesn't seem to be about what is allowed, but what is commonly used. I have a question for you: which gene predictors are YOU using? Or are you using gene annotations produced by others?

Here are a few examples.

NCBI RefSeq

If you're using GFF3 from NCBI, chances are it looks something like this.

NC_007070.3     Gnomon  gene    329235  331223  .       +       .       ID=gene14;Dbxref=BEEBASE:GB42168,GeneID:551678;Name=LOC551678;gbkey=Gene;gene=LOC551678;gene_biotype=protein_coding
NC_007070.3     Gnomon  mRNA    329235  331223  .       +       .       ID=rna27;Parent=gene14;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;Name=XM_624067.4;gbkey=mRNA;gene=LOC551678;model_evidence=Supporting evidence includes similarity to: 66 ESTs%2C 24 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 113 samples with support for all annotated introns;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329235  329459  .       +       .       ID=id117;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    329850  330082  .       +       .       ID=id118;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330166  330301  .       +       .       ID=id119;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  exon    330376  331223  .       +       .       ID=id120;Parent=rna27;Dbxref=GeneID:551678,Genbank:XM_624067.4,BEEBASE:GB42168;gbkey=mRNA;gene=LOC551678;product=receptor expression-enhancing protein 5-like%2C transcript variant X1;transcript_id=XM_624067.4
NC_007070.3     Gnomon  CDS     329333  329459  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     329850  330082  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330166  330301  .       +       0       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1
NC_007070.3     Gnomon  CDS     330376  330416  .       +       2       ID=cds8;Parent=rna27;Dbxref=GeneID:551678,Genbank:XP_624070.1,BEEBASE:GB42168;Name=XP_624070.1;gbkey=CDS;gene=LOC551678;product=receptor expression-enhancing protein 5-like isoform X1;protein_id=XP_624070.1

MAKER

The MAKER annotation workflow (paper, software) is a pretty commonly used gene annotation tool, and produces GFF3 output like this.

scaffold_12     maker   gene    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959
scaffold_12     maker   mRNA    652527  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1;Parent=maker-scaffold_12-augustus-gene-0.959;Name=maker-scaffold_12-augustus-gene-0.959-mRNA-1;_AED=0.24;_eAED=0.18;_QI=0|0|0|0.66|0.5|0.33|3|0|218
scaffold_12     maker   exon    652527  652817  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1203;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    654877  655170  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1204;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   exon    655272  655343  .       +       .       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:exon:1205;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     652527  652817  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     654877  655170  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   CDS     655272  655343  .       +       0       ID=maker-scaffold_12-augustus-gene-0.959-mRNA-1:cds;Parent=maker-scaffold_12-augustus-gene-0.959-mRNA-1
scaffold_12     maker   gene    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851
scaffold_12     maker   mRNA    941547  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1;Parent=snap-scaffold_12-processed-gene-0.851;Name=snap-scaffold_12-processed-gene-0.851-mRNA-1;_AED=0.22;_eAED=0.22;_QI=0|0|0.25|0.25|1|1|4|661|95
scaffold_12     maker   exon    941547  941631  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1206;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942343  942367  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1207;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    942780  942920  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1208;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   exon    943200  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:exon:1209;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     941547  941631  .       +       0       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942343  942367  .       +       2       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     942780  942920  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   CDS     943200  943236  .       +       1       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:cds;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1
scaffold_12     maker   three_prime_UTR 943237  943897  .       +       .       ID=snap-scaffold_12-processed-gene-0.851-mRNA-1:three_prime_utr;Parent=snap-scaffold_12-processed-gene-0.851-mRNA-1

GeneMark, SNAP, Augustus

If we look at the output of several common ab initio gene prediction tools, none of them actually produces GFF3 by default. Here is some output from GeneMark (paper, software).

543   1   +  Initial         799051     799097         47          1 2 - -
543   2   +  Terminal        799236     799266         31          3 3 - -

544   3   -  Terminal        802357     802514        158          3 2 - -
544   2   -  Internal        802607     802685         79          1 1 - -
544   1   -  Initial         802829     802843         15          3 1 - -

Here is some output from SNAP (paper, software).

Einit   8230    8239    +       9.329   0       1       0       scaffold_12-snap.4
Exon    8848    8869    +       2.484   2       2       2       scaffold_12-snap.4
Exon    10121   10208   +       15.302  1       0       2       scaffold_12-snap.4
Exon    11361   11420   +       5.969   0       0       2       scaffold_12-snap.4
Exon    11471   11535   +       0.921   0       2       1       scaffold_12-snap.4
Eterm   12169   12187   +       18.163  1       0       1       scaffold_12-snap.4
Einit   14569   14668   +       -6.918  0       1       0       scaffold_12-snap.5
Exon    15029   15203   +       -1.023  2       2       0       scaffold_12-snap.5
Exon    16171   16307   +       -9.230  1       1       1       scaffold_12-snap.5
Eterm   16667   16698   +       9.829   2       0       0       scaffold_12-snap.5
Einit   17809   17898   -       -5.390  0       0       0       scaffold_12-snap.6
Exon    17299   17350   -       6.978   0       1       2       scaffold_12-snap.6
Eterm   17029   17129   -       -2.857  2       0       0       scaffold_12-snap.6

Both of these formats appear to be tab-delimited, but that's where the similarity with GFF3 ends.

The output of Augustus (paper, software) looks like GTF (a variant of the GFF format) once you lose all the lines beginning with a # symbol, but Augustus also has a GFF3 output mode.

# Predicted genes for sequence number 1 on both strands
# start gene g1
scaffold_12     AUGUSTUS        gene    2841    3400    1       +       .       g1
scaffold_12     AUGUSTUS        transcript      2841    3400    1       +       .       g1.t1
scaffold_12     AUGUSTUS        start_codon     2841    2843    .       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        intron  3027    3097    1       +       .       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     2841    3026    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        CDS     3098    3400    1       +       0       transcript_id "g1.t1"; gene_id "g1";
scaffold_12     AUGUSTUS        stop_codon      3398    3400    .       +       0       transcript_id "g1.t1"; gene_id "g1";
# protein sequence = [MAIKNAEHDLRVIVDAIEGLGLKVAPHKTEAMAFPASALCGRRGAAPPKIRLGGSSILVGSRSRWYISGHSENSSKSP
# RTEGKETTPLQQRDPLDAPLWVSGVVAHCCGGPEGQEGCPGLAAQGSDQGVLRIRDGLLCGYDGCGDHRPRPSDSSAGGGLCRP]
# end gene g1
###
# start gene g2
scaffold_12     AUGUSTUS        gene    4712    15229   0.21    +       .       g2
scaffold_12     AUGUSTUS        transcript      4712    15229   0.21    +       .       g2.t1
scaffold_12     AUGUSTUS        start_codon     4712    4714    .       +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  4858    5591    0.74    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        intron  5686    15028   0.23    +       .       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     4712    4857    0.72    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     5592    5685    0.42    +       1       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        CDS     15029   15229   0.43    +       0       transcript_id "g2.t1"; gene_id "g2";
scaffold_12     AUGUSTUS        stop_codon      15227   15229   .       +       0       transcript_id "g2.t1"; gene_id "g2";
# protein sequence = [MGRNSHRSCCVVNCKITSAKSDCKFYKFPTAKWKINQRKMWVAAVKRQKYIKDEISHAETQTEITEVTGATKVNYANK
# KYICLLFVRTYVLRMLVDVALSNLRFSLFGIRKSLEIFGQSEKADQTRWRLPSCEMEWIESRKGKMRE]
# end gene g2
###

Parsing attributes

So...how should YOU handle attributes in GFF3's 9th column? That depends a lot on what you want to do with the data.

The most important attributes are ID and Parent, which are used to define relationships between features and subfeatures. (These relationships implicitly define a directed acyclic graph of features, although most GFF3 parsers don't directly support traversal of this graph.)

But not only do you have to handle the pre-defined attributes discussed in the GFF3 specification, you also have to be able to handle any number of arbitrary attributes whose keys you may not know beforehand.

The simplest way to handle this would be to parse the attribute column into a dictionary of key/value pairs. Once it's in this form, it's trivial to see what attributes are there and how to access them.

For example, if we have a feature with the following attributes...

ID=mRNA42;Parent=gene19;integrity=0.95;foo=bar

...we would want it in a dictionary like so.

attributes = {
    'ID': 'mrna42',
    'Parent': 'gene19',
    'integrity': '0.95',
    'foo': 'bar',
}

Python code to parse that might look something like this.

attributes = dict()
for keyvaluepair in attributestring.split(';'):
    for key, value in keyvaluepair.split('='):
        attributes[key] = value

One thing to consider is that it's possible for an attribute to have multiple values (separated by commas). It's not commonly used, but it is valid and something that should be handled or at least checked for.


Are there any RepBase alternatives for genome-wide repeat element annotations? I’m using the RepBase libraries in conjunction with RepeatMasker to get genome-wide repeat element annotations, in particular for transposable elements.

This works well enough, and seems to be the de facto standard in the field.

However, there are two issues with the use of RepBase, which is why I (and others) have been looking for alternatives (so far without success):

RepBase isn’t open data. Their academic license agreement includes a clause that explicitly forbids dissemination of data derived from RepBase. It’s unclear to what extent this is binding/enforceable, but it effectively prevents publishing at least some of the data I’m using and generating. This is unacceptable for open science.

Subordinate to this, the subscription model of RepBase also makes it impossible to integrate RepBase into fully automated pipelines, because user interaction is required to subscribe to RepBase, and to provide the login credentials.

RepBase is heavily manually curated. This is both good and bad. Good, because manual curation of sequence data is often the most reliable form of curation. On the flip side, manual curation is inherently biased; and worse, it’s hard to quantify this bias — this is acknowledged by the RepBase maintainers.

Dfam has recently launched a sister resource, Dfam_consensus, whose stated aim is to replace RepBase. From the annoucement:

  Dfam_consensus provides an open framework for the community to store both seed alignments (multiple alignments of instances for a given family) and the corresponding consensus sequence model.

Both RepeatMasker and RepeatModeler have been updated to support Dfam_consensus.

I haven’t tried it yet but it looks promising.


Sequence alignment of extended IUPAC sequences I want to align a short sequence to a reference, and both are encoded in the extended IUPAC notation (i.e. 16 letters).

What alignment software can I use for that?

For example, in Bowtie "all non-A/C/G/T characters in reads become Ns", so it ignores which extended letters are matching.

You can try using Biostrings, depends on how large is your reference, if it is small, you can just pairwiseAlignment with the nucleotideSubstitutionMatrix which supposedly :

library(Biostrings)
data(phiX174Phage)

x = DNAString("TCCAWRSCGCAGAAGT")

lapply(phiX174Phage,function(y)pairwiseAlignment(x,y,
substitutionMatrix = nucleotideSubstitutionMatrix()))

$Genbank
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$RF70s
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$SS78
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$Bull
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$G97
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

$NEB03
Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: -------------TCCAWRSCGCAGAAGT...----------------------------
subject: GAGTTTTATCGCTTCCATGACGCAGAAGT...TAAAAATGATTGGCGTATCCAACCTGCA
score: -21491.92

If it is much larger, you can convert the sequence of interest to a pwm and find it in the reference fasta, like shown in this example


FASTQ file reads extraction paired ends I have paired end reads and want to extract only those which have this sequence within them:
TGTATGTAAACTTCCGACTTCAACTGTA

I tried using grep -B1 -A2 "TGTATGTAAACTTCCGACTTCAACTGTA" file_1.fq | grep -v "^--$"
and the same from the mate pair. After this I tried aligning them , but they dont align because the reads headers are different.
Can someone help in extracting the mate pair of one read from the other file (2.fq), so that they can be aligned ?
regards.

Your question can be rewritten as, "how do I get the union of reads in two files?" Since you have already subset the mates, the steps are:

Process each file, storing the list of reads names in each.
Take the union of these.
Iterate through the original/non-subset files and write only the reads with names in the union to the output.

The files will be both subset correctly and synced properly such that bowtie2 or any other aligner can correctly process them.


How to identify latent variables in single-cell RNA-Seq data I have a single-cell RNASeq sample, in which I'd like to identify latent variables (e.g. response to stress) that I think might be affecting the clustering.

The approach I was planning to use is to reduce the dimensionality of the dataset using one algorithm among:

PCA
ICA
NMF

Then identify the most important biological processes for every component (e.g. in each Principal Component), and then finally regress out the processes that might not be relevant to the analysis or that might affect the downstream clustering.

I have a few questions:

Are there issues in this approach?
Is there a collection of gene sets (similar to MSigDB) that can be used to identify noise latent variables?
Which algorithm (PCA, ICA, NMF, or others) is more suited to this kind of analysis?

I think to identify latent variables, PCA is probably not going to work. NMF might be worth trying. You might want to check out a method called consensus NMF (cNMF) (https://elifesciences.org/articles/43803) published in eLife. The author described some latent variables corresponds to cell types, while others corresponds to cell-cycle or metabolic states of cells.


Generating DNA sequences with constraints I would like some advice on potential strategies to address the following problem.

I want to write a program that will generate DNA sequences that are optimized on two constraints based on an input genome:

The GC content should be as close to the genome as possible.
When taking all k-mers (of a predefined length) in the sequence, the sum of their frequencies in the genome should be as small as possible (i.e. it should be made of rare k-mers).

I currently score sequences using  :
$S = w_{GC} * (1 - |GC_{genome} - GC_{seq}|) + w_{Kmer} * (1 - \frac{mean(f(Kmer_{seq})}{max(f(Kmer_{genome}))}$

Where $w_{GC}$ and $w_{Kmer}$ are user-defined weights defining the relative importance of GC content and k-mer profile) and $f(Kmer_{seq})$ is the frequency of k-mers in the sequence.

Of course I could just use bruteforce and generate many random sequences of appropriate GC content, score them all and pick the very best. But 1) that's no fun and 2) it's slow.

The first thing I attempted was to extend the sequence one basepair at a time. Each iteration has 4 candidates (with a new A, C, T or G added). I pick the candidate using their score ($S$) as probability weight. This "blind" method is too simple and does not yield satisfactory results.

I thought a dynamic approach could work well (filling a table and backtracking),  but I am not sure this is applicable since the cost of extending the sequence would be dependent on all new k-mers generated and not be a fixed value (such as in sequence alignment).

I also thought something like reinforcement learning could be appropriate here, but I don't know much about it, so I am not sure.

I know all of this is probably overkill, but this is mostly for the challenge / learning :)

Which strategy would be most appropriate for this ? Did I miss something simple ? All advice / suggestions are welcome !

You can use a $l$-order Markov chain. Here is the procedure:

Count $l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.
Draw a $l$-mer randomly based on the distribution of $l$-mers. This generates the first $l$ bases.
Let $s$ be the last $(l-1)$ subsequence from the generated sequence. There will be at most four $l$-mers prefixed with $s$. It gives you a distribution of the next A/C/G/T base conditional on $s$. Randomly draw a base from this distribution. You get the next base.
Repeat step 3 until you generate a sequence long enough.

The generated sequence will have the same $l$-mer distribution to the original sequence. The GC content will be the same, too. In your case, if $k$ is small enough, you can let $l=k$. If $k$ is large, you have to choose a smaller $l$.

The procedure above doesn't consider the strand symmetry of DNA, but probably this doesn't matter much. It is also possible to adapt the algorithm to address this.

As to implementation, you can use python, but a high performance language like C will be tens of times faster at least.


Reject reads with low quality bases from a Bam file through pysam I have a code below:

def FilterReads(in_file, out_file):

def read_ok(read):
    """
    read_ok - reject reads with a low quality (<5) base call
    read - a PySam AlignedRead object
    returns: True if the read is ok
    """
    if any([ord(c)-33 < _BASE_QUAL_CUTOFF for c in list(read.qual)]):
        return False
    else:
        return True

_BASE_QUAL_CUTOFF = 30

bam_in = pysam.Samfile(in_file, 'rb')
bam_out = pysam.Samfile(out_file, 'wb', template=bam_in)

for read in bam_in.fetch():
    if read_ok(read):
        bam_out.write(read)

This code works fine by

Rejecting reads having a base with phred quality score below five
But it first takes a BAM file
And then creates a filtered BAM file for further analysis but this takes a lot of time.

So is there a way to reject these reads using pileup in pysam so that I may not have to create a file then sort it and again read it?
Or can I modify this code to perform the same function

For can do this by accessing the basecall qualities from PileupRead.alignment. For example:

_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

This involves checking every read at every position. You could memoise this task:

def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

Or, using functools.lru_cache (thanks @winni2k):

from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff


Why does my SplitDotPlotGG plot show differential expression of a gene but Seurat::FindMarkers() function does not return it? I have analysed  and clustered my single cell rna seq data with methods in the Seurat package. After clustering, I obtained two main clusters: 0 and 1. I used the Seurat SplitDotPlotGG function to visualise the differential expression of a gene of interest X.

SplitDotPlotGG(object,
               grouping.var = "clusterID",
               genes.plot = c("geneX"),
               cols.use = c("blue", "red"), plot.legend = T)

The plot shows clearly that there could be a significant differential expression of geneX between the conditions represented in clusters 0 and 1. However, I cannot see geneX in the list of differentially expressed genes returned when I run:

FindMarkers(object, ident.1 = 0, ident.2 = 1, genes.use = "geneX")

Is there be any explanation for this?
Thanks in advance.

You're showing groups 1 and 2 in the plot but comparing groups 0 and 1 in the differential expression test, so you need to instead set ident.1 = 1, ident.2 = 2 in FindMarkers.

The FindMarkers function by default only returns genes with >0.25 log fold-change expression, with >3 cells in each group. It's possible that these criteria are not being met for gene X. You can try changing these defaults, setting logfc.threshold = 0 to return all genes.


Can I use my computer while MinION is sequencing without negatively affecting the run? I have a Mac Book Pro and I am about to start a sequencing run on the MinION. MinKNOW is going to be running for 2 days.
Can I use my computer during sequencing? Can I browse the web and use Excel, etc? I am not speaking about running a parallelised genome assembly on all my cores using my entire RAM but more mundane activities, like checking my emails or taking notes. Furthermore, can I lock my screen using Control + Shift + OFF without negatively affecting the sequencing?

Yes, that should be fine, especially if basecalling is disabled. Locking a screen should be fine as well.

I prefer to do basecalling after a run anyway, because there's a chance that it could chew up too much processor time and result in lots of skipped reads.


Removing stop codons from alignment file for dn/ds calculation This is one of the question which is very similar to my problem or almost same to my problem

I ran with dummy data nucleotide alignment and the amino file it works perfectly fine,i can run the dn/ds from phangorn without any issue.

But when i ran with my data nucleotide alignment and the amino file I do get this warning

In kaks(carps) : sequence lengths are not a multiple of 3

and error.

The code

Any help or suggestion would be highly appreciated.

I would align your sequences with the software available here: https://github.com/veg/hyphy-analyses/tree/master/codon-msa


Problems replicating Kover paper results I am trying to replicate the results for the creator of this repo:  https://github.com/aldro61/kover2_paper for genotype to phenotype machine learning interpretable (decision tree) prediction.
Specifically in the data section when I try to create the dataset as explained here:
https://github.com/aldro61/kover2_paper/tree/master/data
I make a directory called genomes_dir in the data directory. Then I do:
python download_genomes.py 'mycobacterium tuberculosis/kanamycin/metadata.tsv' genomes_dir

at the end of the download it tells me that it saved file genome_paths.tsv the data directory, not genomes_dir. So I run command:
kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress

and it just hangs with a prompt:
(kover) nobu@gold3forever:~/Desktop/BioInformatics/ResearchAssistant/kover2_paper/data$ kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress
> ^C

it does nothing after hours and hours so I control-C'd to end the process.
Does anyone have any idea what is happening?
I also had to install this to get the datasets from PATRIC:  https://github.com/aldro61/patric_tools

The single quote at 'mycobacterium is not matched so the command prompt is waiting indefinitely for the command to be completed. Maybe you want 'mycobacterium tuberculosis/kanamycin/metadata.tsv' (note closing quote)


How to determine the primary Uniprot accession number from a set of accession numbers? Given a list of Uniprot IDs that are linked to an Ensembl gene ID, is there a way to systematically determine the UniProt accession corresponding to the protein product of the gene's primary transcript with no other information?

According to ExPasy

  Researchers who wish to cite entries in their publications should
  always cite the first accession number. This is commonly referred to
  as the 'primary accession number'. 'Secondary accession numbers' are
  sorted alphanumerically.

But what if the order has been jumbled or compiled from a different source and resorted?

For example:

Ensembl:

ENSMUSG00000035642

Uniprot:

Q8R0P4, Q8CF11, D6RJK8, D6RJJ4, D3Z442, D3Z1Q3, D3YZD8, D3YY39, D3YX09, D3YWY5

This question is cross listed on the Biology stack-exchange.

Reposting my answer from the related ticket in the Biology section:

I think there is an issue with the terminology. The "primary" accession number, is the first accession number in cases where an entry has more than one accession number, as described in http://www.uniprot.org/help/accession_numbers:

  Entries can have more than one accession number. This can be due to
  two distinct mechanisms:

a) When two or more entries are merged, the accession numbers from all entries are kept. The first accession number is referred to as the

  ‘Primary (citable) accession number’, while the others are referred to
  as ‘Secondary accession numbers’. These are listed in alphanumerical
  order.

b) If an existing entry is split into two or more entries (‘demerged’), new ‘primary’ accession numbers are attributed to all

  the split entries while all original accession numbers are retained as
  ‘secondary’ accession numbers.

  Example: P29358 which has been ‘demerged’ into P68250 and P68251.

Both reviewed and unreviewed entries can have primary accession numbers.

What you probably mean, as previous posters understood, are accession numbers of reviewed entries as opposed to unreviewed ones.

In that case, you can indeed add "reviewed:yes" to your query, e.g. when you are using the UniProt ID mapping, http://www.uniprot.org/help/uploadlists


Methods to predict one connected structure from two known separate structures Currently, I am facing such a question: how to predict one connected protein structures from two separate known structures?

It should be different from the protein-protein docking method, which is non-covalent interactions. My question is that the two proteins should "contently" link together.

Any suggestions are appreciated.

There are two approaches here.
Physiological
If domain A and B of your protein interact with each other physiologically, then you need to protein-protein docking. This is not too straightforward as protein are flexible.
Several online servers exists such as Haddock and include ClusPro, which use more traditional approaches. While the online server DEMO (I-Tasser group) does a structural homology based search and a final MD based repacking to better minimise the structures (no protein-protein docking).
Many of these work by magic as far as the user is concerned. But understanding what is going on is really important, even if it is overly laborious. When you don't know where a protein goes, it's called blind or global docking, which is computational expensive and so imprecise that the effort required to prove it is correct for peer-review is too steep for it to be useful. If you are doing local docking it is okay, which may be the case here depending on the linker. Specifically, you need to generate more 'decoys' (potential poses) the more unsure you are: tens of thousands or more. Then you have find clusters and see how each scores following more precise refinement. As an example of a tediously and manually set up run, Rosetta has several docking protocols, which explain the steps required —of note are the many judgement calls, such as in the local refinement, where one has to judge whether setting -dock_pert 3 8 (3 Å and 8° perturbation) is too strict etc.
Pretty picture
The protein-protein docking approach will result in a globular protein, which may be incorrect and majorly is not great for pictures.
If you want a pretty picture of two or more domains most likely you would want them placed like beads on a string. If you have a linker that is not solved the recommended distance to space them is $(3.5 N_{missing})^{0.5}$ Å —unfortunately I do not know where this lore comes from. If it helps I wrote a script that does this in automatic but manual placement works best.
Linkers
In the situation you have a good confidence about how the two parts fit together (or have placed them like beads on a string), you may want to add linkers. If you absolutely don't care about correct bond topology here is a blog post on how to join bits (shoddily) with PyMOL and here is one to do it properly with Rosetta or Pyrosetta.
If one wanted to do a better job, one tool to use is Rosetta remodel, which has a bit of a learning curve, but is very powerful for protein design. This requires the PDB to be on a single chain (pymol alter command, followed by create will do the trick: more here) and actually the amino acids to not have gaps —PDB number being the same as the pose numbering. Rosetta gives better results if the poses are energy minimised, but for this it does not matter too much. Then a "blueprint file" formatted as discussed in the documentation. Briefly, in a way that residues left alone are in the form 1 M . (where residue 1 is a methionine), while the residues before and after the insertion are like 100 D L PIKAA D (meaning D100 is changed to a Asp in loop SS; not actual numbers or residue), while between entries are added 0 X L PIKAA G etc. To model loops, -generic_aa G is a good option as this tells the default rough fitting (centroid mode) to use glycine and not valine as generic residue.
The top pose will actually not have a superb score, so a 3 (regular) or more cycles (15 thorough) with Rosetta Relax application either of the whole protein or with a movemap specific to the linker, may be done.
There actually is an option in Rosetta remodel to have a domain insertion, i.e. it would swing around like a medieval flail. However, it will not repack and therefore properly dock so this option is a bad idea for this specific application.


Merging several fq.gz files or R1 and R2 classes into a single one I have the following files:
V300066187_L4_B5RDBATtnuRAAAAA-408_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-408_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_2.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-407_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-406_1.fq.gz
V300068047_L2_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-407_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-406_1.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_2.fq.gz
V300066187_L4_B5RDBATtnuRAAAAA-405_1.fq.gz

the ending defines whether it is forward _1 or reverse _2.
I want to merge all the _1.fq.gz files into like Merge1_.fq.gz and the same for _2.fq.gz
I tried doing cat *_(number).fq.gz > merged(number).fq.gz and when I used the merges for assembly with MaSurca, it didn't detect any reads and failed. I used one pair and it worked well, but the genome size was smaller than expected, real size = 2 GB, estimated = 1.8GB.
Would appreciate how to solve this issue and use all my short reads for assembly - I need to input one file for all forward reads and one file for all reverse reads.

What you did should work.  You don't need to decompress and recompress.
I'd check the files.  Were they trimmed?  Can you get the top, say 1000 lines of the merged files to work?  (You will need to decompress to count lines accurately)


How to write a hash function for canonical kmers? This question is a follow up from How do kmer counters determine which kmer is 'canonical'?.

In that question we learned that kmer counting programs use a 2-bit hash function to internally represent canonical kmers as they are being counted.

Now I am wondering, how can we implement such a function in C/C++ or in python? More arbitrarily, how could I calculate the canonical kmer hash value using a mathematical function?

For example, how would we transform the 3-mer GAT or the 21-mer GAATACCATAGGATA to 1s and 0s such that:

hash(GAT) == hash(ATC)
hash(GAATACCATAGGATA) == hash(TATCCTATGGTATTC)

You can convert any string hash function to a "canonical" DNA string hash function.

Given a DNA string $s$, let $\overline{s}$ be its Watson-Crick reverse complement. Suppose $h:\Sigma^*\to\mathbb{Z}$ is an arbitrary string hash function. You can define
$$
\tilde{h}(s)\triangleq \min\{h(s),h(\overline{s})\}
$$
Then for any DNA string $s$
$$
\tilde{h}(s)=\tilde{h}(\overline{s})=\tilde{h}({\rm canonical}(s|h))
$$

PS: There are more than one ways to define "canonical" hash function. For example, let $h_0(\cdot)$ be the 2-bit encoding hash function. Given any integer hash function $g(\cdot)$, we can define:
$$
h'(s)\triangleq g\big[\min\{h_0(s),h_0(\overline{s})\}\big]
$$
Then we still have
$$
h'(s)=h'(\overline{s})
$$
You can also replace $\min\{\}$ with XOR or plus to get new definitions.


Comparing phylogenetic models with different datasets I'm a linguist interested in phylogenetic tree inference using language data. I'm posting here because I'm using Bayesian phylogenetic methods in my work (probably using BEAST and/or RevBayes). For the purposes of this question, please accept the assumption that phylogenetic trees are generally applicable to languages (I'm aware this premise may be subject to dispute). My question is intended to be methodological and not necessarily domain specific. You could imagine I'm constructing biological species trees using morphological data.

I have two datasets: A set of binary traits ($D_{b}$) and a set of continuous traits ($D_{c}$).

I'd like to compare the results of phylogenetic tree inference with (a) binary data alone vs. (b) binary data and continuous data together. Model parameters are the same, except of course (b) requires some extra parameters to model the evolution of continuous traits. In other words, both the data and model of (a) is a subset of (b).

I know this probably sounds like a strange research question. The reason it's of interest to me is that the binary data represent a data type that is well-established in historical linguistics and often used to infer language trees. The continuous data come from different language structures and haven't been used to infer language trees before. I want to know if adding this new kind of data can help aid tree inference significantly or if it has a neutral effect (or even just introduces more noise).

I'm familiar with phylogenetic model comparison using Bayes factors. But in the standard use case I'm familiar with, the dataset stays exactly the same and Bayes factors are used to evaluate the fit of different evolutionary models.

My question: Is there a standard way of evaluating the effect of adding more data to a phylogenetic tree model? Can marginal likelihoods be compared still even though the data differs in each model?

If not, would it make any sense to reframe the question as one of model comparison? I'm imagining, for example, running both (a) and (b) with identical datasets ($D_b + D_c$ for both). In (a), there would be a trivial model of trait evolution that doesn't aid tree inference at all—something like a model where traits can jump from any value to any other value with equal probability(?). In (b), there would be a more standard model of continuous trait evolution, e.g. Brownian motion.

The key issue is the data is subject to reticulation, languages are seldom clean breaks. In nucleotide data this can be identified due to deviations from physical linkage and call it recombination. Linkage doesn't apply to languages.

The fundamental assumption is one tree - the same tree - describes all the data, with reticulation that may not be a good assumption. Thus English is a hybrid language, and lots of English words find their way back into French and German 'parents'. Or perhaps a better example is Gaelic mingles with English over time. There are methods used to monitor reticulation in languages.

Thus id keep it simple. The branch length is not only directly comparable it provides linguistical information that is useful to understand. This is the amount of 'evolutionaey' change per word (it that is what you binary data is describing) and the relationship to the node ('breakpoiint') is of primary interest.

The second step is to simply draw a tree of data set A and compare it to data set A+B and data set B and the probabilities describing the sister groups. Are the trees congruent, i.e. are they statistically the same?

The area you describing is more given to parametric bootstrapping, using custom resampling, but you need to be certain that the model being used is appropriate for the data.

The final point is whether a reversible or directional matrix is being used, which impacts on the specific result.


Difference between genome assembly and genome sequence alignment to a reference to find structural variants I'm trying to determine what the difference and benefits of genome assembly and genome sequence alignments are when trying to identify structural variants or transposons in populations.
I've been scouring the internet but have only really come across the difference between short vs long reads and de novo assembly vs reference-based.
My understanding is that to identify variations in structural variants within a population there seem to be 2 main comparative genomic methods, the first being what the 1KGP and SDGP did and sequence the whole genome, align the reads to the reference genome and end up with a BAM file.
The second is to assemble personal genomes and then compare or align the assemblies to each other and the reference genome or using the Lastz/LiftOver/ChainNets Examples: 10.1016/j.gene.2005.09.031
Thanks in advance.

the first being what the 1KGP and SDGP did and sequence the whole
genome, align the reads to the reference genome and end up with a BAM
file.

If you have a well defined reference genome (e.g. human, mouse etc.) and you are interested in population level genetic variation, then this is the main approach. If you sequence a new human genome in the classical way (i.e. short reads ~30X coverage etc.), de novo assembly is pointless generally and you can more rapidly map reads to the reference. The read mapping approach has the large advantage that when you get the to variant calling stage, you can use information about both the depth of sequencing and the base quality scores.

The second is to assemble personal genomes and then compare or align
the assemblies to each other and the reference genome or using the
Lastz/LiftOver/ChainNets

This is the traditional comparative genomics approach used for non-model organisms and for conducting evolutionary rate comparisons at the cross-species level. You discard information on sequencing depth and quality scores, but you can compare many species' genomes at once genome-wide. You also worry less about sequencing errors etc. because you are not looking for relatively rare population level variation, but instead more common (depending on your species) variation between species.


How I swap rownames in these data? I have a microarray gene expression matrix like this by weird gene IDs in rows.

> head(mat[1:10,1:5])

> dim(mat)
[1] 39302    76
>

There is matched gene symbol for each of these identifiers in another matrix

>

> dim(matched)
[1] 23107     1
>

How I can have matched gene symbol with probe identifiers in the row names of my expression matrix please? The problem is, for one gene symbol we may
I tried

> merged <- merge(mat, matched)
Error: cannot allocate vector of size 6.8 Gb

  EDITED

tmp = paste(matched[rownames(mat)],rownames(mat),sep="_")
> rownames(mat) = tmp

EDITED

This is my matrix after matching prob identifiers to gene symbol

> head(array[,1:10,1:5])
               GSM482796 GSM482797 GSM482798 GSM482799
1         OR2T6    0.0171   -0.1100   -0.0394   -0.0141
2          EBF1    0.1890    0.0222    0.0832    0.0459
3 DKFZp686D0972    1.9400    0.2530    0.3770    0.8310
4        ATP8B4   -0.1490    0.0690   -0.0637   -0.0527
5      NOTCH2NL    0.1540   -0.3880    0.2160   -0.0812
6        SPIRE1    0.2920    0.1690    0.5500    0.1430

but now for some genes I have several probes. For example for gene A I have several matched probes. So I have repetition for gene A. How I can take mean over the expression of repeated genes please and having an unique value?

There are actually two questions.

The first one is about the memory issue. I believe merge from data.table would solve that issue.

The second one would be aggregating or summarizing the identifiers. As @StupidWolf mentioned, you need to have a rule, summing them up, averaging them, ... For this second part one way is using dplyr:

merged_data %>%
  group_by(probe_colum) %>%
  summarise_all(mean or sum or ...) %>%
  mutate(new_column = paste(probe_column, gene_column, sep = "_"))

The same thing can be done faster with data.table:

merged_data[, lapply(.SD, mean or sum or ..., na.rm=TRUE), by = probe_column ] [, new_column := paste(probe_column, gene_column, sep = "_"),]

.SD means subset of data and your subsets are your groups based on your probes.


Why I am getting different gene lengths for the same genes in different samples with rsem-calculate-expression? I am running the STAR-RSEM pipeline, I have a simple shell script(*), which I use to loop over fasta files corresponding to different samples. From the output, I am interested in the *.genes.results files, for which the columns are:

> names(my_df1)
[1] "gene_id" "transcript_id(s)" "length" "effective_length" "expected_count" "TPM" "FPKM"

To my surprise, whereas the gene_id and transcript_id(s) columns are the same, I see different length and effective_length values for the same genes in different sample results (omitting the transcript_id(s) column for formatting purposes):

> my_df1 <- fread(paste0("alignment/", my_files[1]))
> head(my_df1[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count   TPM
1: ENSG00000000003.14 2292.23          2242.23            978 26.69
2:  ENSG00000000005.5 1339.00          1289.00             11  0.52
3: ENSG00000000419.12 1024.45           974.45            469 29.45
4: ENSG00000000457.13 4204.00          4154.00            241  3.55
5: ENSG00000000460.16 2953.40          2903.40             62  1.31
>
> my_df2 <- fread(paste0("alignment/", my_files[2]))
> head(my_df2[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count    TPM
1: ENSG00000000003.14 2260.38          2210.38       10052.00 117.78
2:  ENSG00000000005.5 1269.22          1219.22         132.00   2.80
3: ENSG00000000419.12 1058.13          1008.13        4386.00 112.68
4: ENSG00000000457.13 4175.05          4125.05         873.00   5.48
5: ENSG00000000460.16 2783.92          2733.92         780.01   7.39
>
> my_df3 <- fread(paste0("alignment/", my_files[3]))
> head(my_df3[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count  TPM
1: ENSG00000000003.14 2662.59          2612.59             49 2.90
2:  ENSG00000000005.5  940.50           890.50              0 0.00
3: ENSG00000000419.12  956.87           906.87             58 9.89
4: ENSG00000000457.13 6364.00          6314.00             66 1.62
5: ENSG00000000460.16 2661.00          2611.00             15 0.89

The length and effective_length are defined per gene (and its transcripts) and should not differ across samples. Below is from the manual:

  A gene's 'length' and 'effective_length' are defined as the weighted
  average of its transcripts' lengths and effective lengths (weighted by
  'IsoPct')

I would appreciate any hints pointing out what I am missing.

* The script I am using:

for file in ...; \
do rsem-calculate-expression \
-p 8 \
--star \
# --output-genome-bam \
--star-path path-to-star-executable \
--star-gzipped-read-file \
--forward-prob=0 \
$file \
path-to-genome \ # should end with the string used for the "reference_name" parameter of "rsem-prepare-reference"
${file:11:7}; \ # sample prefix
done

EDIT: For those who will be interested, this entry from the RSEM Users Google Group demonstrates how to calculate gene length from weighted transcript length.

If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.


How can I find the relevant pathway map from gene-gene or protein-protein interaction list? I have some difficulties to understand/interpret the pathway map and how a gene-gene interaction list or DNA sequencing can map into pathways.

In addition what's the difference between MARK/ERK pathway and the output of KEGG

See a project like Plant reactome, where they infer the pathways of all sequenced plant genomes using orthology to well annotated species.
http://plantreactome.gramene.org/


Searching for gene expression data by cell line I have two cancer cell lines (OCI-Ly18 & riva) that I want to find gene expression data for, but I'm not aware of many gene expression databases that allow searching by cell-line without searching by gene.

I tried Genevestigator on the recommendation of this thread
but found no data on either cell line.

What would be the best way to find gene expression data for these cell lines?

Try the Gene Expression Omnibus - it looks like they have some datasets.


E.coli Sequencing & Analysis I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.

The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.

I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.

For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.

Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!

Thanks,

Alex.

A de novo assembly with spades should work as well as a mapping with bowtie2.
While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.
When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.
I'd like to give you the advice to use some kind of distribution system like Anaconda which helps you to keep your system clean and make it easier to reproduce your steps.
Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :
De novo assembly
1. Quality Check ( you already did this )
2. Run Spades
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log

This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like.
You now have a folder called k12 which holds the assembly information.
To get some information about the contigs you can use the emboss toolkit.
infoseq k12/contigs.fasta

fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 -              NODE_859_length_78_cov_1.000000 -              N    78     43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 -              NODE_860_length_78_cov_1.000000 -              N    78     44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 -              NODE_861_length_78_cov_1.000000 -              N    78     43.59

This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.
3. Remove low quality contigs
There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.
4. Check for possible contamination
So to check for possible contamination you can use kraken2. But it needs some disk space. I am using here the minikraken database
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log

The ecoli_kraken.report shows an overview of how the contigs could be classified to availible species from the database.
Based on this results you can also remove some contigs.
5. Run Prokka
Now after you removed unwanted contigs, you can run prokka to get the annotation.
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log

And you will find the information in the prokka folder
Mapping
The mapping process is a little bit easier. You can use bowtie2 and samtools
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam

158216 reads; of these:
  158216 (100.00%) were paired; of these:
    35696 (22.56%) aligned concordantly 0 times
    120047 (75.88%) aligned concordantly exactly 1 time
    2473 (1.56%) aligned concordantly >1 times
    ----
    35696 pairs aligned concordantly 0 times; of these:
      4458 (12.49%) aligned discordantly 1 time
    ----
    31238 pairs aligned 0 times concordantly or discordantly; of these:
      62476 mates make up the pairs; of these:
        56512 (90.45%) aligned 0 times
        5365 (8.59%) aligned exactly 1 time
        599 (0.96%) aligned >1 times
82.14% overall alignment rate

You now also get some information about the similarity to your reference sequence.
Now you need to work a little bit more on the BAM file. You need to sort and index this.
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam

Now you can create a consensus sequence using bcftools and seqtk
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta

This gives you now a sequence where unknown bases are marked as N.
As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file bio.yaml
name: bio
channels:
  - bioconda
  - defaults
dependencies:
  - aragorn=1.2.38=h470a237_2
  - barrnap=0.9=2
  - bcftools=1.9=h4da6232_0
  - bedtools=2.27.1=he941832_2
  - blast=2.6.0=boost1.60_0
  - bowtie2=2.3.4.3=py36h2d50403_0
  - emboss=6.6.0=h6debe1e_0
  - hmmer=3.2.1=hfc679d8_0
  - infernal=1.1.2=h470a237_1
  - kraken2=2.0.7_beta=pl526h2d50403_0
  - libdeflate=1.0=h470a237_0
  - libidn=7.45.0=2
  - minced=0.3.2=0
  - parallel=20160622=1
  - perl-app-cpanminus=1.7044=pl526_1
  - perl-bioperl=1.6.924=4
  - perl-carp=1.38=pl526_1
  - perl-constant=1.33=pl526_1
  - perl-data-dumper=2.161=pl526_2
  - perl-encode=2.88=pl526_1
  - perl-exporter=5.72=pl526_1
  - perl-extutils-makemaker=7.34=pl526_2
  - perl-file-path=2.15=pl526_0
  - perl-file-temp=0.2304=pl526_2
  - perl-parent=0.236=pl526_1
  - perl-threaded=5.22.0=13
  - perl-xml-namespacesupport=1.12=pl526_0
  - perl-xml-parser=2.44=pl526h3a4f0e9_6
  - perl-xml-sax=1.00=pl526_0
  - perl-xml-sax-base=1.09=pl526_0
  - perl-xml-sax-expat=0.51=pl526_2
  - perl-xml-simple=2.25=pl526_0
  - perl-yaml=1.27=pl526_0
  - prodigal=2.6.3=1
  - prokka=1.12=3
  - samtools=1.9=h8ee4bcc_1
  - spades=3.12.0=1
  - tbl2asn=25.3=0
  - boost=1.60.0=py36_0
  - bzip2=1.0.6=h14c3975_5
  - ca-certificates=2018.03.07=0
  - certifi=2018.10.15=py36_0
  - curl=7.62.0=hbc83047_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - giflib=5.1.4=h14c3975_1
  - icu=58.2=h9c2bf20_1
  - jpeg=9b=h024ee3a_2
  - libcurl=7.62.0=h20c2e04_0
  - libedit=3.1.20170329=h6b74fdf_2
  - libffi=3.2.1=hd88cf55_4
  - libgcc=7.2.0=h69d50b8_2
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgd=2.2.5=hceca4fd_3
  - libpng=1.6.35=hbc83047_0
  - libssh2=1.8.0=h1ba5d50_4
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libtiff=4.0.9=he85c1e1_2
  - libuuid=1.0.3=h1bed415_2
  - libwebp=1.0.0=h222930b_1
  - libxml2=2.9.8=h26e45fe_1
  - ncurses=6.1=hf484d3e_0
  - openjdk=8.0.152=h46b5887_1
  - openssl=1.1.1=h7b6447c_0
  - perl=5.26.2=h14c3975_0
  - pip=18.1=py36_0
  - python=3.6.7=h0371630_0
  - readline=7.0=h7b6447c_5
  - setuptools=40.5.0=py36_0
  - sqlite=3.25.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - wheel=0.32.2=py36_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=ha838bed_2

Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.
conda env create -f bio.yaml

I hope this gives you some help how to process your data. If you need some information just let me know ;)


Calling for one specific SNP from multiple sequencing runs Just learning the basics of bioinformatics and bash shell scripting. Is there a correct way to search for a specific allele with one or even multiple sequencing runs? I am searching for prevalence of the C2416T allele in the SARS-CoV-2 genome as a mini project. I downloaded a bunch of sequencing runs, I think I need to convert them to bam file type than I am not sure how to call for one specific SNP, especially in a loop for 5-7 runs. Can this be done with bash?
I was looking into variant calling from a bam file type using vcf and I can't find a way to call for one specific snp on multiple bam files at once. Ideally I would want to filter the variant out of the .bcf files or use grep since I sort of understand those the most and I’m not working with a large amount of data or a human genome.

Getting SNPs from whole-genome data is usually done via mapping sequencing reads on the reference (bwa-mem or bowtie are two popular mappers) and then using one of variant calling tools (like freebayes, samtools, or GATK). However, I would also look up if there is something virus-specific.
According to this answer, you could get variant calls of a specific region using samtools if you give it a bed file that contains the region of interest (three columns chr    start    stop format). For example using human reference
samtools mpileup -uDl regions.bed -f hg19.fa file1.sorted.bam file2.sorted.bam | bcftools view -bvcg - > RAL_samtools.raw.bcf

You could also specify a bunch of files, not just two and it would do the variant calling on them all.


A reliable fetcher of short read using SRA/ENA accession I am trying to build a workflow that gets data automatically from databases of sequencing reads (SRA - sequence read archive, ENA - european nucleotide archive). Till now I was pulling everything from SRA and my script have worked, but now I encountered a first entry that is only in ENA (ERR2135453), not in SRA. And I found that my download script has its limitations.

My bash script looks like this (requires aspera ascp and fastq-dump) and I wanted to make changes that are commented now

#!/bin/bash

ACCESION=$1
OPENSSH=/my/part/to/key/asperaweb_id_dsa.openssh

# if [ ${ACCESION:0:3} == SRA ]; then
    URL=anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/
# else
#    URL=ena_url
# fi

ascp -QT -i $OPENSSH $URL

fastq-dump data/"$ACCESION".sra --outdir data/ --split-files --gzip
rm data/"$ACCESION".sra

The problem is in getting ENA url. Unlike SRA url that I can easily get from the accession (anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/${ACCESION::3}/${ACCESION::6}/"$ACCESION"/"$ACCESION".sra data/), the ftp server of ENA is bit more confusing. The url corresponding to ERR2135453 accession is ftp://ftp.sra.ebi.ac.uk//vol1/err/ERR213/003/ERR2135453. Where 003 there comes from? Also not all reads are in err subdirectory.

Is there an easy way how to download reads from ENA knowing just the accession number?

--- edit ---

I found on the relevant ENA site about downloading data using aspera that ftp urls should be possible to derive only from accession also in ENA.

  Submitted read data files are organised by submission accession number
  under vol1/ directory in ftp.sra.ebi.ac.uk:

  ftp://ftp.sra.ebi.ac.uk/vol1//

  where  contains the first 6 letters and
  numbers of the SRA Submission accession.

  For example, the files submitted in the SRA Submission ERA007448 are
  available at: ftp://ftp.sra.ebi.ac.uk/vol1/ERA007/ERA007448/.

The reason why there is this strange 003 will probably be that it's quite recent submission and it has not been properly categorized yet (other accessions I tried were on "expected" urls).

However, the problem still persist. Is there a way how to fetch data of new submissions using accessions only?

Apparently there is no 100% reliable solution, but in the most cases, you should be able to use the enaBrowserTools python scripts (enaDataGet) to quite reliably download data.

The run ERR2135453 is a special case because the sample submission (ERS1939590) for it has been cancelled. The solution for these reads is to download the data manually from the ftp location.

This answer is based on kind assistance of Suran from ebi helpdesk.


Is it possible for coronavirus or SARS to be synthetic? I have heard several conspiracy theories regarding the origin of the new coronavirus, 2019-nCov. For example that the virus and/or SARS were produced in a laboratory or were some variant of Middle Eastern respiratory syndrome (MERS), shipped via laboratory workers.

I am well aware bioinformatics has debunked many conspiracy theories involving infectious diseases, an important one being polio vaccination programs in Africa were the origin of HIV, for example here. Likewise, numerous conspiracy theories on the man-made origins of HIV were similarly overturned using bioinformatic studies.

Do you have bioinformatics evidence that would debunk the current conspiracists about coronaviruses?

The scenarios are impossible and would be laughable if they were not so serious. The evidence is in the phylogenetic trees. Its a bit like a crime scene when the forensics team investigate. We've done enough crime-scenes often going to the site, collecting the pathogen, sequencing and then analysis - (usually neglected diseases) without any associated conspiracy theories.

The key technical issue is coronaviruses are zoonoses , pathogens spread to humans from animal reservoirs and phylogenetic tree really helps understand the  how the virus is transmitted.

Trees

The key thing about all the trees are bats. Bats lineages are present at every single point of the betacoronavirus phylogeny (tree), both as paraphyletic and monophyletic lineages, one example is this tree of betacoronaviruses here. Meaning the nodes connecting the branches of the tree to the "master-branch", represent common ancestors and these were almost certainly bat-borne cornaviruses. This is especially true for SARS and - here bat-viruses are EVERYWHERE.
The tree here also shows that SARS arose on independently on two occassions, again surrounded by bat lineages and 2019-nCov has emerged separately at least once, again associated with bats.
Finally the tree below is a figure from BioRxiv Zhou et al (2020) "Discovery of a novel coronavirus associated with the recent pneumonia outbreak in humans and its potential bat origin" shows the 2019-nCov lineage is a direct descendent of a very closely virus isolated from a bat (RaTG13*). This is a really conclusive finding BTW.

Note, I don't normally present inline images, but it is such a nice finding (hint to reviewers) and BioRxiv is open access.

Conspiracy theory 1: laboratory made virus

Literally it would require someone passaging a new virus, with unknown human pathogenicity, and independently introducing all the earlier passages enmass across bat populations of China. They would then hope each lineage becomes an indepedent virus population before, then introducing the virus to humans. Thus when field teams of scientists go around using mist nets to trap the bats, buy them from markets, isolate the virus and sequence it they would find a beautiful, array of natural variation in the bat populations leading up to the human epidemics, that perfectly matches vast numbers of other viral zoonoses. Moreover, this would have to have happen substancially prior to SARS and 2019-nCov, because the bat betacoronaviruses have been known about prior both epidemics, viz. its simply not feasible.

Biological explanation
General Bats are a reservoir host to vast numbers of pathogens, particularly viruses, including many alphaviruses, flaviviruses, rabies virus and beleived to be important in ebolavirus (I don't know about this) and even important to several eukaryotic parasites. It makes sense, they are mammals, so evolutionary much closer to us than birds for example, with large dispersal potential and roost in 'overcrowded' areas enable rapid transmission between bats.

Technical
The trees show bats are the common ancestor of betacoronaviruses in particular for the lineage leading into the emergence of 2019-nCov and SARS, this is seen in this tree, this one and the tree above. The obvious explanation is the virus circulates endemically in bats and has jumped into humans. For SARS the intermediate host, or possible "vector" was civet cats.

The theory and the observations fit into a seamless biological answer.

Conspiracy theory 2: Middle Eastern connection

I heard a very weird conspiracy theory attempting to connect MERS with 2019-nCov. The theory was elaborate and I don't think it is productive to describe here.

Biological explanation
All the trees of betacoronaviruses show MERS was one of the earliest viruses to diverge and is very distant from 2019-nCov, to the extent the theory is completely implausible. The homology between these viruses is 50%, so its either MERS or 2019-nCov. Its more extreme than mixing up yellow fever virus (mortality 40-80%) with West Nile virus (mortality <<0.1%), the two viruses are completely different at every level.

What about errors? Phylogeneticists can spot it a mile off. There are tell-tale phylogenetic signatures we pick up, but also we do this to assess 'rare' genetic phenomina. There is nothing 'rare' about the coronaviruses. The only anomaly is variation in the poly-A tail and that is the natural variation from in vitro time-series experiments. Basically we've looked at enough virses/parasites through trees, that have no conspiracy theories at all (often neglected diseases), and understand how natural variation operates - so a phylogenecist can shift the wheat from the chaff without really thinking about it.

Opinion
The conspiracy theories are deeply misplaced, and the only connection I can imagine is its China. However, the Chinese have loads of viruses, influenza in particular which causes major pandemics, but that is a consequence of their natural ecology (small-holder farming) allowing the virus to move between reservoir hosts. I've not visited small-hold farms in China, but I have in other parts of the world and when you see them, you get it. The pigs, chickens (ducks .. China), dogs, horses and humans all living within 10 meters of each other.

Conclusion
Shipping large numbers of bats to market, bat soup, raw meat from arboreal (tree-living) mammals such as civets that are sympatric to bats. Then consider the classical epidemiology in light of the phylogenetic data, which is very consistent, a single picture emerges that coronavirus is one of many zoonoses which has managed to transmit between patients.

Summary The fundamental point is the bioinformatics fit into the classical epidemiology of a zoonose.

*, Note The bat coronavirus RaTG13 predates the 2019-nCov outbreak by 7 years. It is not even clear whether the virus has been isolated, i.e. could just be a RNA sequence.

"They have found some 500 novel coronaviruses, about 50 of which fall relatively close to the SARS virus on the family tree, including RaTG13—it was fished out of a bat fecal sample they collected in 2013 from a cave in Moglang in Yunnan province."

Cohen, "Mining coronavirus genomes for clues to the outbreak’s origins" Feb, 2020 Science Magazine,


Which measure should be used in a PCA or RNA-seq data? TPM or counts? I'm trying to understand the magnitude of batch effects in my RNA-seq samples, and I was wondering which expression units are more suitable to draw a PCA. I'm thinking of either counts or TPM, but things like rlog or vst could work too.

Additionally, I'm wondering whether either of this units should be log-transformed first, to avoid high-abundance transcripts driving the PCA.

tldr: log transform counts and TPMs, but rlog/vst are preferred

TPM should be log transformed to get more useful results. If you're using DESeq2 already (given the reference to rlog and vst, this seems likely), then please go ahead and use rlog or vst. That will give you more reasonable results than raw counts. If you're stuck with counts for some reason, then first use normalized counts so they're at least a bit more comparable and then log transform them so your highly expressed genes aren't driving everything.

Edit: As an aside, if you know what the batch effect is (e.g., library prep. date), it's sometimes convenient to include that in your model. You can then assess the genes that are actually changed due to that, which is sometimes handy to know (e.g., which genes might be more/less prone to degradation).


Can I ask DESeq2 which variable alone explains which gene's behavior the best? Imagine I have a dataset of highly periodic gene expression, taken at densely spaced timepoints. I know that the expression is periodic, but different genes have different offsets.

It is trivial to convert my timestamps into a cosine (cos_time = cos(offset + time*pi*2/PERIOD)), so that this new variable is linear with the genes' expressions -- but only with the ones that have a peak at offset.

What if I could have multiple cos_times (let's call them phases) and then ask DESeq2: "OK, which gene's behavior is most linear with which phase?" Barring comparing p-values of each LRT test, is there a better way?

Disclaimer: I know about cosinor, Lomb-Scargle periodograms, ARS, JTK etc, and I know they are likely more applicable in this particular case. I'm interested in whether it is in general possible to ask such a question from DESeq2, even if the variables in question are not periodic. Phase just seemed like the simplest concept to demonstrate this on.

If you have densely spaced time points the using non-linear regression to directly fit the signal to a cosine or sine curve, which would directly produce the offset, would be much simpler. DESeq2 wouldn't be the appropriate tool then, you'd use base R and the nls() function.


Read sorted and indexed BAM files faster in C++? I have some sorted and indexed alignment large BAM files. I'd like to read all the reads, and apply custom operation on each of them. The orders aren't important. I'm using htslib in C++ to do the reading. I have a single powerful machine.

The problem is ... reading a large BAM file is very slow even in C++.

Q: Is there any trick to speed up the reading? Reading BAM files multithreaded (I have many cores)?

In this answer I want to show you some benchmarks that compares three different serial ways of reading data in C++, the third one being the fastest.

In the first method I use an std::istreambuf_iterator, in the second I read the file line by line into a std::vector and in the third one I use C flavoured operations.

/*
compile with g++ -std=c++11 -O2 reading_files.cpp -o reading_files
test the program with the following
dd if=/dev/random of=./dummy5GiB bs=1024 count=$[1024*1024*5]
 for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
*/

#include <fstream>
#include <iostream>
#include <vector>
#include <cstring> //for strcmp
#include <string>
#include <algorithm> //for min

using namespace std;

int main(int argc, char * argv[])
{
    if(argc != 3) {
        cerr << "usage: reading_files <path> <number 0 or 1 or 2>\n";
        return 1;
    }

    ifstream file(argv[1], ios::in|ios::binary);
    if(!file.is_open()) {
        cerr << argv[1] << " file not opened!\n";
        return 1;
    }
    cout << "---------------------------\n";
    if(!strcmp(argv[2],"0")) {
        cout << "method 0\n";
        vector<char> content((std::istreambuf_iterator<char>(file)), std::istreambuf_iterator<char>());
        cout << content.size() << " bytes read\n";
    } else if(!strcmp(argv[2],"1")) {
        cout << "method 1\n";
        string line;
        vector<string> content;
        while(getline(file,line)) {
            content.push_back(line);
        }
        cout << content.size() << " lines read\n";
    } else if(!strcmp(argv[2],"2")) {
        cout << "method 2\n";
        file.seekg(0, std::ios::end);
        unsigned long long file_size = file.tellg();
        cout << "file size is " << file_size << "\n";
        //maximum 5 GiB of memory, so to fit my RAM
        unsigned long long max_block_size = 1024ULL * 1024ULL * 1024ULL * 5;
        file.clear();
        file.seekg(0, std::ios::beg);

        unsigned long long block_size = min(max_block_size,file_size);
        char * buffer = (char*)malloc(block_size*sizeof(char));
        unsigned long long processed = 0;

        while(!file.eof()) {
            file.read(buffer, block_size);
            cout << "read a block of " << min(file_size - processed, block_size) << " bytes\n";
            //do work here with the block in memory (if the block size is > 0)
            processed += block_size;
        }

        free(buffer);
    }
    file.close();
    return 0;
}

In the third method I read the file in blocks of 5GiB, so to fit the free RAM of my machine.

Here is the benchmark of the program on a 5 GiB randomly generated file.

bash-3.2$ for i in 0 1 2; do time ./reading_files dummy5GiB $i; done
---------------------------
method 0
5368709120 bytes read

real    0m35.194s
user    0m19.665s
sys 0m12.999s
---------------------------
method 1
20973135 lines read

real    0m40.478s
user    0m34.046s
sys 0m5.870s
---------------------------
method 2
file size is 5368709120
read a block of 5368709120 bytes
read a block of 0 bytes

real    0m4.413s
user    0m1.757s
sys 0m2.617s

Here is the benchmark on a real 32GiB .bam file using only the third method, here the splitting of the file read plays a fundamental role. I can't tell you how big because my hard disk memory gets saturated before the end of the execution and so I have to kill the program.

time ./reading_files ./HG00252.mapped.ILLUMINA.bwa.GBR.low_coverage.20130415.bam 2
---------------------------
method 2
file size is 34316054058
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 5368709120 bytes
read a block of 2103799338 bytes

real    0m50.656s
user    0m9.135s
sys 0m35.728s

The key point is to minimise the file operations reading large block of data at once (as long as they fit the memory) and operate on them.

Using mmap could lead to an improvement, but I don't think it would be as large as using the third method instead of the first two.


What format is this? Pretty sure it's not a BED file Time for guess the bioinformatics file format... :)

I found this '.bed' file on GEO (GSE84660), but it's clearly not a BED file. Anyone know what it might be? And what can view it? It's from a HiC expt.

#column headers: bait1_chr, bait1_start, bait1_end, bait2_chr, bait2_start, bait2_end, contact_ID, NA, NA, B1_D0_counts, B1_D3_counts, B1_D6_counts, B2_D0_counts, B2_D3_counts, B2_D6_counts

chr1    831895  848168  chr1    850619  874081  b2b_1   .       .       114     125     127     75      97      92
chr1    831895  848168  chr1    889424  903640  b2b_2   .       .       15      12      16      4       15      20
chr1    831895  848168  chr1    903641  927394  b2b_3   .       .       13      13      10      18      15      8
...

It looks like a weird way to represent data. You're right, it does not at all look like a BED.

Digging in GEO made me find this information:

  Supplementary_files_format_and_content: Tab-delimited file reports one
  contact per row. Interacting HindIII fragments are represented in bed
  format (chr/start/stop) with bait (or upstream bait, in the case of
  b2b contacts) listed first. The following fields correspond to a
  unique contact ID, fragment IDs for first and second fragments, and
  raw read counts supporting contacts in each replicates.

So I think the authors merged a BED3 format plus additional BED-unrelated fields generating a tab-separated file - not clearly stated in their methods though. The easiest would be to contact the authors maybe?

Depending on your aims, if you only need the first 3 fields (but I doubt it) I would reformat this "bed" file (to make it processable as a BED file) in two sub-files:

cat original.bed | cut -f 1,2,3 > bait1.txt
cat bait1.txt | head -3
chr1    831895  848168
chr1    831895  848168
chr1    831895  848168

cat original.bed | cut -f 4,5,6 > bait2.txt
cat bait2.txt | head -3
chr1    850619  874081
chr1    889424  903640
chr1    903641  927394

Not sure if it would fit your aims though.

Also, for visualization of contact matrices, maybe you'll find this paper useful.


How to compact variant data to their genes without bias? I have a dataset of genes I am trying to collect data on from public databases, to use as features in machine learning. I am trying to take some features from UCSC genome browser (e.g. number of CpG islands per gene, number of DNase clusters per gene, regulatory enrichment scores etc.) however I am not sure how to control for bias where a gene that is larger in length - and so will then have more CpG islands or higher regulatory enrichment scores simply due to gene length.
Is there a way to correct for gene length when taking/condensing variant data to individual genes?
Edit:
For reference, my machine learning model aims to predict whether a gene is the most likely to be causal for a disease (out of all the genes given to the model). The model will score the genes as a regression classification between 0 to 1 (0 being least likely to cause disease and 1 being most likely to cause disease). I plan to later further investigate the genes with the highest scores.
The model uses a variety of multi-omic features (e.g. GTEx gene expression the genes have for many tissues, GWAScatalog data, gene intolerance scores, protein-protein interaction data, drug interaction data, phenotypic scores etc.). However, I am missing epigenetic data to describe my genes so I've been looking to collect based on UCSC's variant data (CpG islands, histone modifications, DNase clusters) - however this leads to my gene length problem when I am trying to reliably take data from the variant level.
I've been plotting my features and gene length, and seen that the UCSC epigenetic data does correlate with having a larger gene length if there is a higher count of regulatory sites (0.8 r2 for some), and so this is what I'm looking to correct.

Its very easy, just let the ML sort this out for you and that is its advantage, You're thinking of GLM style calculation where you pre-screen the data with bivariate plots, where there needs to be nice Q-Q plots and low residual.
For ML simply include the gene length as one of your parameters along with CpG etc ... and the ML regression analysis SVC, lasso, ridge, random forest will figure the relationship out between gene length and CpG. You do zero, the ML does everything, hence from a statistical point of view purists object because you don't know the relationshiop the ML has deduced between the variables, but you will get regression weights for non-DNN stuff, which will give you some idea of the impact of length.
There is the issue of transformations and that can be complicated, but I'd try untransformed data first. The only disadvantage of this approach is the user will have to input the gene size when they want to check out your training algorithm.


Is the number of contigs and scaffolds an artifact of genome assembly and rather arbitrary, or is it canonical for each species? E.g. bowhead whale genome has X contigs and Y scaffolds, but couldn't these easily be different numbers?

This is purely arbitrary and an artifact of how well the assembly turned out. Better quality assemblies with have fewer and longer contigs.


TopHat2 versus HISAT2 inner workings In my intro to bioinformatics course, we mentioned that TopHat2 and HISAT2 will both try to align as many reads as possible to the reference genome (TopHat2 has been superseded by HISAT2). For the reads that were not mapped, these probably cross exon boundaries. To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments. Then:
"TopHat2 tries to align many fragments from each read using all possible exons as the search space. HISAT2 focuses on one fragment from each read, and once that fragment has been anchored it only searches locally."
My question is about that last part "once that fragment has been anchored, [HISAT2] only searches locally." Isn't it possible that a given fragment will map to many different regions of the genome? So how is it possible that it will be "anchored"? Are we going to anchor it multiple times and see which anchor makes most sense or what exactly?

I don't believe this statement is true:

To solve this issue, both TopHat2 and HISAT2 chop up the unmapped reads into around 25bp fragments.

TopHat2 did not support spliced alignment directly.  This is why it followed this general approach of chopping the read up into small parts, followed by "contiguous" alignment of the sub-reads.  Hisat2 / Hisat-genotype support spliced alignment directly insofar as they don't pre-partition the reads into small parts, but allow dynamic extensions of alignments across splice junctions (in a way similar to how STAR works).  This is generally a fundamentally better approach and is one reason, among many, why TopHat/TopHat2 have been put into maintenance-only mode and have been officially superseded by Hisat2/genotype.


Generating the reconstructed alignment from BAM I have a (small) BAM file with CIGAR and MD fields.

Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...

My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.)

Here is read and reference:

CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
|||||||||||||||| |||||||||||||||||
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG

Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :

ReadAligned    Reference
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
...    ....

Question 2:
Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library cigar or the cigar parsing functions of pysam

If all you want are the number of insertions/deletions, you can use this to extract that info:

from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

Deletions would be the same thing with 'D' in place of 'I'


How to make chromosome color maps for bed ranges I have genomic .bed file data of 4 different types; type A,B,C,D. These are some examples-

Type A:

1   101380000   101710000   A
1   110085000   110320000   A

Type B:

1   100930000   101335000   B
1   10430000    10560000    B

Type C:

1   101775000   101825000   C
1   103570000   104070000   C

Type D:

1   153325000   153475000   D
1   156355000   156415000   D

This data is for chr1, but I such data for all the chromosomes. I want to draw a chromosome map (possibly using the Gviz package?) which can assign different colors to the 4 types of genomic ranges and make me an idiogram like plot for the chromosomes I wish to see, so I know where my ranges physically are, on the chromosome.

Can this be done using the Gviz package in R? If not, what is the best way and code to do this?

You can use karyoploteR to plot your regions on an ideogram quite easily.

Disclosure: I'm one of the authors of the tool.

For this example I'll start creating 4 GRanges with random regions, but with your data you should load your bed into R and split it into four GRanges.

library(karyoploteR)

typeA <- createRandomRegions(nregions = 100, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 20, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 1000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 100, length.mean = 2e6, length.sd = 5e6)

colors <- c(A="#FFD700", B="#00BFFF", C="#FF7F50", D="#90EE90")

And we can plot them. We have to start creating an empty ideogram with plotKaryotype (we are working with human here, but you can use any genome you want) and then plot the regions in the GRanges with kpPlotRegions. As a bonus, if two regions overlap it will plot them one above the other.

kp <- plotKaryotype(genome = "hg19")

kpPlotRegions(kp, data=typeA, r0=0, r1=0.2, col=colors["A"])
kpPlotRegions(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"])
kpPlotRegions(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"])
kpPlotRegions(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"])
legend(x = "right", legend = names(colors), fill = colors, border=NA)

And you'll get something similar to this

In this example we have fairly large regions and not a lot of them. As @heathobrien pointed out, plotting the actual regions in the whole genome plot might not be the best visualization option. You can either use zoom to create different plots of smaller parts of the genome or change your visualization to show the density of regions along the genome with kpPlotDensity.

Plotting the density and tweaking a bit the visualization (and changing the example to have 10 times more regions) with this codei

typeA <- createRandomRegions(nregions = 1000, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 200, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 10000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 1000, length.mean = 2e6, length.sd = 5e6)

kp <- plotKaryotype(genome = "hg19", plot.type = 4,
                    ideogram.plotter = NULL, labels.plotter = NULL, main="Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)

kpPlotDensity(kp, data=typeA, r0=0, r1=0.2, col=colors["A"], window.size     = 10e6)
kpAddLabels(kp, labels = "Type A", r0=0, r1=0.2)
kpPlotDensity(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"], window.size = 10e6)
kpAddLabels(kp, labels = "Type B", r0=0.25, r1=0.45)
kpPlotDensity(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"], window.size = 10e6)
kpAddLabels(kp, labels = "Type C", r0=0.5, r1=0.7)
kpPlotDensity(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"], window.size = 10e6)
kpAddLabels(kp, labels = "Type D", r0=0.75, r1=0.95)

it would look something like this

The plots could be tweaked a lot changing various parameters and you could add other data on top of it if needed. You can find more information on the Bioconductor landing page or at the karyoploteR's Tutorial and Examples page.


Is Wikidata a reliable KB for genomic and drugs annotations? I'm looking for a Linked Open Data approach to annotate a dataset with human genome and drugs information.

According to the Linked Open Data cloud, there are a lot of interconnected RDF vocabularies devoted to life sciences but, unfortunately, a lot of them seem out of date.

I recently found Wikidata:

  Wikidata is a free and open knowledge base that can be read and edited
  by both humans and machines.

It seems exactly what Linked Data should be and it seems continuously updated. At a first sight, it covers very well genomic and drugs knowledge.

Does Anyone know if Wikidata is a reliable source for genomic and drugs annotation?

The question is what you mean with "reliable". For bacterial genome annotations there is one Wikidata based platform called WikiGenomes (publications).


nanopore - where to retrieve information from the basecaller used Where can I find information about the basecaller used in a specific nanopore run (we have a Gridion machine)?

I know already that we have Guppy v. 3.0.3, which I found out by clicking on the about tab on the MinKNOW GUI).

I was wondering, though, if additional information about the basecaller is outputted in any of the files generated when the run is finished? I checked already the final_summary.txt and it does not contain such information. I also checked the head of the fastq files and no further details.

My ultimate goal for the moment is to run medaka on a draft genome assembly. But for that I'd need to specify the medaka model. From medaka_consensus -h I see that the options available and defaults are:

-m  medaka model, (default: r941_min_high_g330).
        Available: r941_min_fast_g303, r941_min_high_g303, r941_min_high_g330, r941_prom_fast_g303, r941_prom_high_g303, r941_prom_high_g330, r10_min_high_g303, r10_min_high_g340, r941_prom_diploid_snp, r941_min_high_g340_rle.
        Alternatively a .hdf file from 'medaka train'.

Therefore I already know that my option should contain min and g303 (format of the models available described here). But know I'd need to find out the {pore} and {caller variant}

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" .


How to understand methylation level? For Bisulfite Sequencing data, we can calculate methylation level for a specific locus,
$$
\mbox{methylation level} = \frac{\mbox{# methelyated reads}}{\mbox{# total reads}}
$$
Here I think, for this locus, it can only have one status, methelyated or not. Or to say, methylation level here can only be 0 or 1. Why does it have other values in (0,1)?

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.


How to represent a deletion at position 1 in a VCF file? I am writing a small script to write a VCF file from short multiple sequence alignments, and stumbled into an interesting corner case.

This is the way a deletion should be represented in VCF format (I assume version 4+):

NC_000016.9  2138199 rs1799758   GTGAG   G  .   .  RS=1799758;...

But what happens if the deletion is at the first position? Immagine the following two aligned sequences:

>reference
AGTC
>sample
-GTC

How should it be encoded in VCF format? According to this thread on SEQAnswers, it might be necessary to include the first position after the end of the deletion, like this:

chromosome  1  AG G  .  .  [...]

Is this a VCF 4+ correct way to encode such variants?

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


Receptor-Ligand R-package for visualization I was wondering if there was any R-package for receptor-ligand interaction similar to SingleCellSignal (https://github.com/SCA-IRCM/SingleCellSignalR_v1/blob/master/SingleCellSignalR/vignettes/UsersGuide.Rmd)

+++ if there are vignettes and/or the graphs can be customized.

+++ if it can be linked to Seurat

Note: I have tried iTALK, but I keep getting an error that seems to be a common problem and have left a query.

Thanks

Two nice tools for receptor-ligand interaction analysis are:

CellPhoneDB by the Teichmann lab, although this is not an R package...
NicheNet by the Saeys group

Note that these tools are mainly meant for analysis of receptor-ligand interactions, although I think NicheNet also has some visualization options. CellPhoneDB needs to be run from a shell, but the results can easily be imported into R for visualization. You might check the papers on these 2 methods for some ideas of visualization. Combining Seurat's visualization options and ggplot2 should get you already pretty far.