How to plot genomic.fna fasta file in R using Gviz?

I am not familiar with Gviz, but I am familiar with a variety of tools for visualizing genome features. At the scale of entire microbial genomes, the primary sequence itself is an overwhelming amount of information and can't really be visualized in a meaningful way by itself. I mean, you could color the nucleotides, but few patterns in DNA can be discerned by an unassisted human eye.

What we're typically interested in are features encoded in the genome sequence: protein-coding genes, non-coding RNAs, transposable elements, regulatory elements, binding sites, epigenetic modifications, and so on. A particular scientist is usually only interested in a subset of these features, based on what particular questions they are investigating. In your comment you say you "would like to annotate protein coding regions and compare them to other genomes", so it seems like we're on the same page here.

Visualization tools often display the genome linearly (even, for example, if it's a circular bacterial chromosome), zoomed out so that the nucleotide sequence isn't visible, and overlayed with boxes, arrows, or other glyphs representing the genes or other features of interest. Here are some examples:

This question thread contains images generated with the GenomeTools AnnotationSketch program (gt sketch command).
The UCSC genome browser which provides web-based access to dynamic displays of the human genome and many other mammalian genomes.
The IGV browser can be downloaded and easily run on your personal laptop or desktop computer.
I'm not familiar with Gviz, but it definitely looks like it belongs in this class of tools.

In each of these cases, the visualization software is not doing the work to determine the positions of the genes and other genome features. This was done independently using tools designed for this task, storing the information in an annotation file in GFF3, GTF, or BED format. The visualization tools need both the genome sequence (in Fasta format) and one or more annotation files (in GFF3, GTF, or BED format) to make meaningful visualizations.


How to install/import Biopython into Python 3.8/ PyCharm IDE

Install Anaconda3 the instructions are here. Then making sure you are connected to the internet type
conda install -c anaconda biopython

Once Anaconda is installed don't attempt to use pip, if you do then you need to set an 'environment' in Anaconda and work from that environment.
I suspect you used pip or pip3 via
pip install biopython


How to extract exome on-target reads in batch?

It sounds like bedtools intersect will work for you:

bedtools intersect -wa -a <alignment.bam> -b <region.bed> > <intersect_alignment.bam>

  When a BAM file is used for the A file, the alignment is retained if overlaps exist, and exlcuded if an overlap cannot be found.  If multiple overlaps exist, they are not reported, as we are only testing for one or more overlaps.

It's easy enough to wrap this up into a for loop for batch processing (line breaks can be removed if desired):

 for alnBam in alignmentFiles*.bam;
   do echo ${alnBam};
   bedtools intersect -wa -a ${alnBam} -b <region.bed> > intersect_${alnBam};
 done

There are a lot of other options which should cover most related uses. For more details, just run "bedtools intersect" (with no arguments).


about two group comparison of three group data in DESeq2 package

As explained on the comments, having some other samples might distort your estimations, add an outlier or something alike.

The advantage of having all the three (or all) groups of samples before testing differences, because usually methods like limma, DESeq2 and others uses the information of variation of the genes across all samples to update their prior probabilities and the variation of the genes on the statistical model used to compare between groups.

If the method and your model are robust use all the three groups (generally the methods used are robust, but sometimes the models are not). By robust I mean that it doesn't miss any know (and available) source of variation. Otherwise the comparison might be comparing only the samples not the effect of the condition on the samples. If you only compare the samples it would be equivalent to use a t.test to the limma methods (not completely sure on the DESeq2 case). Or explained differently, you wouldn't take advantage of the statistical improvements implemented on these methods

Usually the reasons why these groups are sequenced on the same batch is because they are related (or to add some known ground truth), so you can release it with the data when communicating your results but you may omit that you used this third group on the methods section.


Difference between PPI networks and Gene Co-expression networks

Based on my rudimentary understanding, PPIs are constructed from gene co-expression networks

Not really, no. PPI stands for protein-protein interaction. In PPI networks, each node is a protein and each edge (line) connecting it to another node represents an interaction. A very common example is that an edge represents a physical interaction between these two proteins. This means that there have been experiments to indicate that these two proteins come into contact in the cell and interact.

For example, they might form part of a protein complex, so their interaction was seen using chromatin immunoprecipitation, where the two proteins were captured together. Or, they may have been seen together using FRET, a technique that lets you detect when two molecules are in very close proximity.

Alternatively, some interactions can be theoretical, where in silico modelling has suggested that these two proteins can interact, even though we don't have any direct evidence that they do so.

Finally, yes, it is also possible to make a PPI network based on gene co-expression data, under the assumption that if two genes are co-expressed, then their protein products will be present in the same tissue. However, that doesn't mean they interact in any way, so it isn't really a very common way of building PPI networks, as far as I know.

In co-expression networks, on the other hand, each node is a gene and not its protein product. And here, the edges indicate the coexpression. Often the level of correlation between the two genes' expression data.

The most important difference between these two kinds of network, apart from the obvious mentioned above, is that they are used in different ways. A gene coexpression network lets you identify genes that are regulated in similar ways. It is a tool you can use to investigate the regulatory network of the cell, to see what genes are "turned on" or "off" together.

The PPI network, on the other hand, is a way of understanding the result of this gene co-regulation. You are not looking at how the gene network is regulated any more, but how the products of the genes interact with one another. So it is a very different focus.

You can maybe think of this in terms of school children and their parents. You could make a network where each node represents a pair of parents, and then draw lines between them if their children were born in the same year. That would give you a general overview of when this group of people procreated. However, that two children were born on the same day doesn't necessarily imply that the children themselves will ever come into contact. This could be the equivalent of the gene co-expression network.

Then, make another network where each node is a child and draw a line between any children who know each other. This would be the equivalent of the PPI network. You would likely find some correlations between the two, with children who know each other being more likely to have been born at around the same time, but the two networks give you very different information.


Detection of CNV(InDel) of intermediate size

My impression is that small InDel (a couple of bp) is identified through cigar string in BAM...

Well, yes, by gaps in the alignment, easily detected by examining the CIGAR string.

  ...and typical CNV (at least thousands of bp) is detected through read depth.

Also by deviations from expected read pair distance or orientation.

  What about InDel or CNV with size these them, say hundreds of bp? These CNV is too long to be covered by read and too short to be detected by statistically comparing the difference of read depth. Do we have consensus that this kind of CNV is hard to detect?

Yes, what you describe has certainly been recognized and acknowledged in the literature. It was one of the primary motivations for a recently published paper I contributed to on de novo (germline) variant discovery. Here's an excerpt from that paper:

  In a reference-mapping context, calling indels with confidence requires accurate mapping of each read spanning the indel, with all gaps arranged consistently. This is possible only for short indels and tends to be prone to error and misalignment. Thus prediction of indels with length 10 bp has proved to be very challenging and accompanied by high false-positive and false-negative rates. Furthermore, the prediction of SVs via read mapping is only possible through indirect signatures such as alterations in read depth or read-pair signatures. These signatures can be quite noisy and result in high rate of false-negative and false-positive prediction.

As for your question:

  is there a common strategy of detection for them?

Numerous mapping-free k-mer based methods for variant detection are coming up in the literature, and all seem to perform pretty well. The trick is how to define which k-mers are interesting—then you grab those k-mers (or the reads containing them) and assemble them into variant-spanning contigs. A few examples:

In the Kevlar paper, we defined  "interesting" k-mers as those with high abundance in an individual of interest but zero (or very low) abundance in the individual's parents.
The HAWK paper, which describes a k-mer based GWAS-like method, defines "significant" k-mers as those significantly over- or under-represented in one class of samples versus a second class.
The NovoBreak paper compares paired tumor and normal WGS data to find k-mers that span mutatino breakpoints in cancer genomes.

  Do we generally think them having impact on protein function (especially in oncology background)?

There is definitely potential that indels can have a substantial impact on protein function. But since their detection has been so difficult for so long, a lot remains unknown about their actual collective impact on function.


How can I calculate the silent mutation for each position in codon?

Calculate in this context would just mean determine. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.

BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. Terming things like this as "silent polymorphism" or "silent change" or something along those lines is preferable.


Does MinKNOW work with Mac OSX high sierra 10.13.1?

Some replies from ONT is:

  OS X is supported and you can view the installation instructions on
  our Downloads page linked below. At the moment High Sierra is not
  officially supported as we have not yet fully tested it, but other
  users in the community have reported no issues with it so far.

  The compatibility test only performs a simplistic check to see if you
  are using one of our officially supported OS X versions. Since you are
  on High Sierra, it will show a failure. That does not mean the
  software can't be installed and used. As I said, other users have
  reported no problems, but we can't officially support it if you
  encounter an issue.

So, I am going ahead and installing MinKNOW on my High Sierra.


Repetition-free Longest Common Subsequence as similarity measure

A brief literature review can give you the answers.

This preprint gives the definition of RFLCS more clearly:

  A repetition free Longest Common Subsequence (LCS) of two sequences x and y is an LCS of x and y where each symbol may appear at most once.

Adi et al proposed to solve the genome rearrangement problem raised in this paper. If you want to know more, follow the references and citations of these papers and read them.

Basically, Adi et al think the length of RFLCS measures the distance between two genomes. The longer the more similar. I am not sure if this is a good measurement, though.


API/website for blasting peptide sequence against database of all plant proteomes

A simple BLAST on the UniProt website does the trick (target database plants).

Best matches are subsequences ETVEREKE and EKKRRETV in Green alga and Eelgrass respectively.


How to convert mzXML to mzML format (mass spectrometry files)?

ProteoWizard offers this functionality:

msconvert your_mzXML_file --mzML

FileConverter of OpenMS (I believe using the integrated ProteoWizard) also does this.


annotation using ChIPseeker package error


Help with the definitions of fields in a VCF output by Strelka

What is the difference of ReadCount in INFO column and DP in FORMAT column?

I'm not super familiar with Strelka but I believe it follows the VCF spec, so DP in the FORMAT column should give the depth of coverage (ie read count) per sample. ReadCount in INFO is not in the VCF spec but based on the definition it sounds like the depth of coverage aggregated over all samples.

  What is the difference of number of allele versus number of reads covering the variant position ?

I believe the "number of alleles used" refers to counts of each base that are used in the indel error modeling process as described in the preprint here:

https://www.biorxiv.org/content/10.1101/192872v2.full

So it sounds like subsets of the number of reads covering the variant position.


merging two/or more bed file into one bed file

In my opinion, the easiest way to merge bed files is to use bedtools merge.


Are sequencing reads written differently when inserts and indexes are sequenced seperately?

The order of sequencing in a dual index flow cell is:

Read 1
Index 1 (I7)
Index 2 (I5)
Read 2

Between each step there is a switch of the primer used (Read1-primer -> Index1-primer (I7) -> Index2-primer (I5) -> Read2-primer)
There are differences in how the I5 is dealt with between dual- and single-read flow cells as explained in the Illumina Indexed Sequencing Overview. In this document there is a lot more detailed explanation of the sequencing workflow that may be helpful to you.
But as explained by the other answers, none of this will matter in most cases. Even if you decide to output the Index sequences they will be exported as separate FASTQ files.


How should I address batch effects in my experiment?

Yes, though thankfully Your 2018 PreD samples will help you resolve this.
Simply add the batch effect to the design (~Batch + Treatment) and DESeq2 (or edgeR or Limma) will handle this for you.

You do not need SVA or RUV, thankfully, since you quite cleverly sequenced one group in both batches.

To clarify, your coldata will be something like:

Group    Time    Batch
A        Pre     A
A        Pre     A
A        Pre     A
A        Post    B
A        Post    B
A        Post    B
...
D        Pre     A
D        Pre     A
D        Pre     A
D        Pre     B
D        Pre     B
D        Pre     B

And the design is actually ~Batch + Group + Time. You can probably add an interaction term in if you then modify the model matrix manually to ensure it's full rank.


What RNA-Seq expression value would be closest to Microarray equivalent?

I think it is very hard to say which are the closest because they are not really comparable. But since you are using Spearman correlation, I guess RPKM, FPKM, and TPM do not change the order of gene expression levels. You might also want to normalize RNA-seq and microarray data so that they are more comparable.


Keep Format and Individual fields when annotating VCF with VEP

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


Single sample in group: normal pipeline or Kal's Z test

You'll have a bit more power if you use all of the samples at once.
Nothing is useful for this, don't waste your time with it.
You're not looking for outliers, you're looking to see how samples cluster. You can just use correlation for that.


Standard Way to Preprocess Gene Expression?

First off - do not use FPKM for across-sample normalization. This is such a common issue that I have a prepared text for it:

Why not to use FPKM

They have a poor sensitivity when looking at low expressed genes [1, 2]. Short transcripts/genes may have falsely over-estimated expression values [5].
They may not normalise adequately across samples [3]. The reason for this is: FPKM are normalised across libraries using a single variable, namely the “total number of reads”. This library size factor estimation can easily be biased by highly expressed, highly variable genes, unlike the method used by the other approaches, which model the library size factor based on summary statistics extracted from all genes. In other words as stated in [4]: “The median ratio method is usually quite robust. The biggest problem for library size correction is not actually the genes with low counts, but the genes with very high counts, which have high variance and disproportionately influence estimators like the total sum. This is why the total sum is a notoriously bad estimator for library size.” Or the same by Lior’s Pachter himself [6].
They have a high False Positive Rate in differential expression analysis [1, 2]. A lot of the genes/transcripts deemed to be differentially expressed are not, mostly due to the effects of point 1), i.e. high artefactual variability in the expression and point 2) i.e. inadequate library size estimation.

On the other hand, approaches such as DESeq, DESeq2, vst or voom have adequate between-sample normalisation, have a controlled FDR and can reliably determine differential expression even for lowly expressed genes if there is enough replication or if an adequate log2 fold change cutoff is selected [7]. If you really insist on using cufflinks, consider performing upper quartile normalisation, which corrects for the point 1) issue.

References
[1] Dillies M.A., Briefings in Bioinf. 2012
[2] Soneson & Delorenzi, BMC Bioinf. 2013
[3] Lior Pachter’s lecture at CSHL (Cold Spring Harbour Laboratory) from minutes 32 on: https://www.youtube.com/watch?v=5NiFibnbE8o
[4] From Answer: "Over-correction" in the size-factors of the DESeq2 package:
[5] http://www.biomedcentral.com/1471-2105/14/370
[6] https://liorpachter.wordpress.com/2014/04/30/estimating-number-of-transcripts-from-rna-seq-measurements-and-why-i-believe-in-paywall/
[7] http://rnajournal.cshlp.org/content/22/6/839

What I recommend, is that you get the raw count data from TCGA, create a count matrix using the DESeq2 package in R, use the function varianceStabilizingTransformation() from DESeq2 to further process your counts and model the networks from there. Your pipeline already mentions voom, which is another good alternative, as is limma. All of these are available on Bioconductor.


Distinguishing Splicings that Produce Functional Transcripts Amongst Alternates Shown

There are two resources worth checking: GTEx portal and Ensembl.
Uniprot and NCBI will call "canonical" the longest sequence comprising valid exons. These are not necessarily the biologically valid ones.

According to Ensembl (TP53, there are a bunch that are experimentally validated and are all class 1. Note the Uniprot column being nearly homologeous: to match the transcripts one has to look at the length, in this case TP53-201 is 393aa long, which is class 1. As mentioned this is not always the case.

So which are more expressed in tissue X?
In GTEx, under the collapsed div labelled "Exon Expression" press the "isoform Expression" tab. The most expressed is ENST00000620739.4 (354aa long) join with ENST00000269305.8 (393aa).

In the Exon Expression you will find out what exons are the issue.

Now. to get the sequence of a transcript or covert a variant Ensembl drops the ball. There is no clear web tab that I know of. Nor an API route.
ftp://ftp.ensembl.org/pub/release-99/fasta/homo_sapiens/pep/Homo_sapiens.GRCh38.pep.all.fa.gz contains all the sequences however. For a project of mine (ENSTMapper in this link) to map an amino acid mutation on one transcript to another I ended up write my own short converter, although there are a few similar tools out there.


Nextflow: ordering the execution of processes; outputs of process A are inputs of process B

Nextflow uses channels (asynchronous FIFO queues) to communicate between processes. The interaction between these processes, and ultimately the pipeline execution flow itself, is implicitly defined by these input and output declarations.
The issue with your code, I think, is the use of multiple input channels. Consider the following example that uses only a single queue channel. To avoid declaring the inputs again in process A's output declaration, we can output a key to join on as the first element in the output tuple. In the case below, this is simply the 'some_name' variable. We can then join up the two channels using the join operator as follows:
Channel
    .fromFilePairs( './data/*.txt', size:1 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .into { input_ch_A; input_ch_B }

process A {

    input:
    tuple val(some_name), path(some_path) from input_ch_A

    output:
    tuple val(some_name), val("procA_${some_name}"), path("procA_${some_name}") \
        into procA_results

    """
    touch "procA_${some_name}"
    """
}

process B {

    input:
    tuple val(some_name), path(some_path), val(some_other_val), path(some_other_path) \
        from input_ch_B.join(procA_results)

    """
    true
    """
}
```


Although generally I recommend people to be wary about the warnings that FastQC gives (it tends to be overly paranoid), the GC content graphs here do look odd. It's good that you've had a look at it with another program to confirm the observation. They should have a distribution that is fairly close to normal, a bit like this:

Without more information about your sample, I can't tell whether or not this is expected from your input. This looks like fairly substantial contamination, but could also be caused by biased selection during sample preparation (e.g. low quality sample, low yield, PCR amplification bias), or just looking at a multi-organism sample rather than a single-organism sample.

As a contamination check, I'd recommend running your RNASeq reads through a kmer-based classification tool (e.g. OneCodex, Centrifuge, BBSketch, or Sourmash). This should tell you the most likely place that contaminant sequences came from.


The UniProt website offers an ID mapping tool that tends to work pretty well:

Then, hit "submit" and you will get the list of corresponding PDB IDs:


Verify a predicted protein in one genome in a different genome of the same species

This should be very easy to do. Here are some options:

Use a tool like Exonerate or GeneWise that can align protein sequences to genomic DNA while attempting to model splice sites etc.
As you said, blast. There's nothing special in what you describe, really, this is what tBLASTn is for. Just use your protein as the query sequence and tBLASTn as the blast flavor which will take a protein query and match it against all 6 possible reading frames of a DNA database. Set up your assembly as the blast DB, and then blast as you would normally.

Of the two, using the first approach is much better since it will try to build a valid gene model for you instead of simply finding the regions of high homology as blast will do.

Finally, just for completeness, you can go old school and run a full de-novo gene prediction tool on your assembly. Something like GeneID for instance. But that's certainly overkill if you're just looking for a single gene.


Phylogenetics analysis between diverse species

I second commenter Devon Ryan's suggestion, I'm just going to give a little more detail about how to do that. You can plug the BED files into the UCSC browser following these instructions: https://genome.ucsc.edu/goldenPath/help/customTrack.html

Now you can look at the overlap between your bed tracks (black boxes I believe) and conservation (red box).

Assuming that is all you want to do, that should do it. If you want to do further data analysis or processing, you can download the conservation track as a .wib file from UCSC. You can convert this file to BED as described here, and use BEDTools to analyze it together with your BED of regions of interest.


How to know which ortholog database is the best for my work?

To recover orthologous sequence for model species
You should be able to find your orthologs using BioMart and your gene IDs.
There are numerous way to do it, directly on the ensembl website or through API. You can start from these links (quite old though) :

http://www.ensembl.info/2009/01/21/how-to-get-all-the-orthologous-genes-between-two-species/
https://www.youtube.com/watch?v=EwfOQsr2O-M #Get orthologs via BioMart
https://www.biostars.org/p/12245/ #Finding Human/Yeast Orthologs Using Embl?

Which database should I use
It is broad and depends on the nature of your work. In your case, working with "major model species", DIOPT (the first line of the table in Quest For Orthology) has the very feature you need on its front page for going from mouse to human. More eukaryotic/viral oriented database will be used for taxonomic identification, annotation inference, etc.
I hope it helps.


Set directory where MinKNOW writes FAST5 files

I found a page explaining where this is set up in a user_conf file here:

https://community.nanoporetech.com/support/faq/test1/minknow/minknow/how-do-i-change-the-directory-for-the-reads-folder?search_term=user_conf

Windows

Open Notepad with administrative rights (search for notepad, right click and then select run as administrator). Select open file and navigate to C:\ProgramFiles\OxfordNanopore\MinKNOW\conf and then open user_conf to view the script. Under output_dirs change the two instances of "C:\data" (default read location) to the directory required. Please note the "\" is important.

Mac OS

Go to /Applications/MinKNOW.app/Contents/Resources/conf/user_conf and open user_conf to view the script. Under output_dirs replace the default path "/Library/MinKNOW/data" with your path of choice.

Ubuntu

Go to /opt/ONT/MinKNOW/conf and open user_conf to view the script. Under output_dirs replace the default path "/var/lib/MinKNOW/data/reads" with your path of choice.


What is a standard approach to binarize microarray gene expression data?

Based on the info you provide, ArrayBin R package provides you the necessary tools:

binarize.array() from ArrayBin, allowing:

  Implementation of an adaptive method for binarizing gene expression data on a per-probe basis and demonstrate the superior effectiveness of our method when compared with other, commonly used approaches. This adaptive binarization method can be applied to DNA methylation microarray data,
  which has implications for cross-platform integration, and can reduce batch effects in the data.

In complement, I'd advise you go through some literature and choose the approach which fits you best. Some suggestions:

Computational analysis of microarray data
Binary matrix factorization for analyzing gene expression data
and probably most relevant for you given your application: Multiscale Binarization of Gene Expression Data for Reconstructing Boolean Networks

Maybe you also want to look into more general concepts of transformation prior to microarray data analysis:

Microarray data normalization and transformation
A Tutorial Review of Microarray Data Analysis


That's a good question. I remember seeing a poster about making trees directly from reads, but unfortunately, I can not find any other record than this tweet. Perhaps you can try to email David, if he would be keen to share the implementation.

While searching for any information about David's work, I run into a GitHub repository called reads2trees. It's quite old and it's not super popular, so hard to say if it can be useful. But the name matches the problem description quite well.

I also remember there were some attempts (An example paper) to get phylogeny directly from kmers, but I have never tried that on my own. I also think these approached were focused on microbial phylogienies. So again, not sure if it's going to be useful for you.

That's all I got.


How do I get the compartment Id and volume from SBML using libSBML Python?

There aren't any examples in the python docs on how to do this. The best bet is to search the C/C++ docs which are more extensive. There is a python example in the docs which shows how to retrieve the number of compartments in a model but it doesn't go further. The code below shows how you can retrieve the information requested, assuming you have a model called testmodel.xml:

import libsbml

document = libsbml.readSBML("testmodel.xml")

if document.getNumErrors() > 0:
   print("Encountered the following SBML errors:" + "\n")
   document.printErrors()
else:
   model = document.getModel()

   nCompartments = model.getNumCompartments()

   print ("Number of compartments: ", nCompartments)
   for i in range (nCompartments):
       comp = model.getCompartment(i)

       print ("Id = ", comp.getId())
       print ("Volume = ", comp.getVolume())


Convert Reactome Protein IDs to Pathway IDs?

If you don't mind hitting it 50k times and are OK with python3...

from urllib import request
import json

def getPathways(proteinID):
    baseURL = 'http://reactome.org/ContentService/data/query'
    PathwayIDs = set()
    try:
        response = request.urlopen('{}/{}'.format(baseURL, proteinID)).read().decode()
        data = json.loads(response)
        if 'consumedByEvent' in data:
            for event in data['consumedByEvent']:
                PathwayIDs.add(event['stId'])
        if 'producedByEvent' in data:
            for event in data['producedByEvent']:
                PathwayIDs.add(event['stId'])
    except:
        pass
    return PathwayIDs

Usage would then be something like:

l = ['R-HSA-49155', 'R-HSA-199420', '']
for rid in l:
    ids = getPathways(rid)
    for _ in ids:
        print("{}\t{}".format(rid, _))

Which would produce:

R-HSA-49155 R-HSA-110239
R-HSA-49155 R-HSA-110240
R-HSA-49155 R-HSA-110238
R-HSA-49155 R-HSA-110356
R-HSA-199420    R-HSA-8948800
R-HSA-199420    R-HSA-6807106
R-HSA-199420    R-HSA-6807206
R-HSA-199420    R-HSA-6807126
R-HSA-199420    R-HSA-8847968
R-HSA-199420    R-HSA-8850997
R-HSA-199420    R-HSA-2321904
R-HSA-199420    R-HSA-8948775
R-HSA-199420    R-HSA-8944497
R-HSA-199420    R-HSA-6807134
R-HSA-199420    R-HSA-8850945
R-HSA-199420    R-HSA-8873946

Note that this will silently ignore invalid or missing IDs such as '' (that's the try and except above. Note also that these are different pathway IDs than what you provided in your example. The main reason is that the protein IDs you showed are not always involved in the pathways IDs you showed (in my example, they always are).


Missing data mappings in mygene.info while trying to convert Genes Ensembl Ids to Entrez Ids

As others have mentioned in their answer, bioMart is usually the best place to go for this infomation because it draws its data directly from the Ensembl database.

However, you will find you do not get a full translation of Ensembl IDs to Entrez IDs even there. The reason for this is there simply does not exist a simple 1-to-1 mapping of Entrez IDs to Ensembl IDs: they are not just different names for the same thing, but are actually different things. At some point someone (or more likely some algo) has decided that the structure Ensembl ID A is referring to is similar enough to the structure that Enztrez ID B is referring to to link the two together, but such a linking is not always possible.

In the example you gave, the Ensembl database does not list ENSG00000267809 as being equivalent to 4730, neither does the Entrez database list the opposite, even though both genes have the same name. Because of this, it is very unlikely that any service you use to convert the IDs will pick up this match.

In this particular case, the reason for this is probably that the gene is a processed pseduo-gene and is not annotated in Refseq/NCBI as producing any transcripts.


What format is this? Pretty sure it's not a BED file

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.


Structural variant calling for low-coverage PacBio data

There is an evaluation of PB Honey and Sniffles algorithms for low coverage PacBio datasets in this preprint and another evaluation is shown on this poster. Both reports agree that optimal is (surprisingly) combination of PB Honey and Sniffles.

Author of Sniffles have benchmarked Sniffles against PB Honey, where he shown that Sniffles performs significantly better. Take a look on this presentation (slide 15).

Another option is SMRT-SV, but I am not aware of any benchmarking.


Efficiently aligning a lot of reads on the same small reference sequence

I would suggest using a whole genome aligner like BWA or BOWTIE2 and swapping out the whole genome reference with your specific reference.

These alignment algorithms are designed to align large numbers of NGS reads against a single reference.  They can be used with any reference including short ones like your example or small viruses.


What do you call screening results?
Do you have a kind of array with protein or RNA data?
If so, you can use several open-source tools available on the internet, like Enrichr or network analyst to better understand the pathways affected by your treatment.


how to train a gene dataset with a nearest shrunken centroid classifier?

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

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

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

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

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


Beside the correct question by @Kamil, this sounds like a job for datamash.

First let's get the number of columns in the file

$ head -n1 genotypes.txt| awk -v FS="\t" '{print NF}'

Now lets calculate the average per column with datamash

$ datamash -H mean 2-<n> < genotypes.txt > mean.txt

Replace <n> with the number of columns we received with command above.

For your example the output will look like this:

mean(SNP1)  mean(SNP2)  mean(SNP3)  mean(SNP4)
0,5 0   1,5 1,5

One can remove the mean with sed

$ sed -i -e 's/mean(//g' -e 's/)//g' mean.txt


What are the advantages and disadvantages between using KEGG or Reactome?

One big downside of KEGG is the licensing issue. One big advantage of Reacome are various crosslinks to other databases and data.

ad 1, This depends on which pathway, they are both primary databases. Sometimes other databases that for instance combine data of primary databases have better annotation of pathways (there is an example in the review paper bellow)

ad 3, There is very extensive relatively new (2015) review on this topic focused on human pathways: Comparison of human cell signaling pathway databases—evolution, drawbacks and challenges. However I could not find there which one is more complete ...


Loading all PPI information from Uniprot into a graph in NetworkX

The query link for all interactors would be thus : https://www.uniprot.org/uniprot/?query=organism%3A9606+AND+interactor%3A*&sort=score

There is always an option of downloading the static version of the service if you find it incovnient to deal with their API for your use case from here and have it all to yourself, pipe data into your dataframes locally. I find it useful to just broswe their repositories sometimes instead of reading the specifications.

In general, you can get a variety of datashapes  from uniprot that you cannot get from their static database-image by specifying the appropriate parameters in the request. You can browse and explore parameters by expanding their GUI's search fields as described in this little guide. Then you can request your data like so, in python for example.

import requests
URL = "https://www.uniprot.org/uniprot/"
PARAMS = {'organism': 9606, 'format': 'tab', ....#Additional parameters#...}
response = requests.get(url=URL, params=PARAMS)


Modifying protein_scale in ProteinAnalysis in SeqUtil in Biopython?

There is not a way to control the edge weights as written, but making a modified version is straightforward. Namely make a class that inherits the class Bio.SeqUtils.ProtParam.ProteinAnalysis and overwrites the method relevant methods.
Having a gander at the files one sees the following:
def protein_scale(self, param_dict, window, edge_list: list):
    """[...]
    """
    # generate the weights
    #   _weight_list returns only one tail. If the list should be
    #   [0.4,0.7,1.0,0.7,0.4] what you actually get from _weights_list
    #   is [0.4,0.7]. The correct calculation is done in the loop.
    weights = self._weight_list(window, edge)
    scores = []
    ...

So the weights are generated by the 'private' method ._weight_list. Underscore means it is for internal use and for production code it is bad to assume that the name of the private method many not change, but Biopython has not changed in years so you'd be safe. And this is for a quick hack anyway...
class JankyProteinAnalysis(Bio.SeqUtils.ProtParam.ProteinAnalysis):
    overriding_weights = None

    def _weight_list(self, window, edge):
          if self.overriding_weights is not None: # dodgy override!
               print("Weight method hack done!")
               return overriding_weights
          else: # behave normally
               return super()._weight_list(window, edge)

So you'd do:
Y = JankyProteinAnalysis('MSTAEFAQLL')
Y.overriding_weights = [0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.]
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,float('nan'))

Note that the edge argument is no longer doing anything.


Use conditional in workflow in Nextflow dsl2

Is this inside the main.nf file or inside the process.nf file? Anyway, the param.anything variable should already be available in the nextflow context when executed, thus, you could simply do:
if (params.aligner == "bowtie") {
    align_bowtie2(get_fq_info.out.fq_info)
} else {
    align_bwa(get_fq_info.out.fq_info)
}

Check out this main.nf of a nextflow pipeline of mine, it contains some working if/else snippets.

Here it is another StackOverflow question on conditional nextflow pipelines that might help.

Hope it helps.


Say you are sequencing to 2X coverage. Suppose at a site, sample S has one reference base and one alternate base. It is hard to tell if this is a sequencing error or a heterozygote. Now suppose you have 1000 other samples, all at 2X read depth. One of them has two ALT bases; 10 of them have one REF and one ALT. It is usually improbable that all these samples have the same sequencing error. Then you can assert sample S has a het. Multi-sample calling helps to increase the sensitivity of not so rare SNPs. Note that what matters here is the assumption of error independency. Ancestry only has a tiny indirect effect.

Multi-sample calling penalizes very rare SNPs, in particular singletons. When you care about variants only, this is for good. Naively combining single-sample calls yields a higher error rate. Multi-sample calling also helps variant filtering at a later stage. For example, for a sample sequenced to 30X coverage, you would not know if a site at 45X depth is caused by a potential CNV/mismapping or by statistical fluctuation. When you see 1000 30X samples at 45X depth, you can easily know you are looking at a CNV/systematic mismapping. Multiple samples enhance most statistical signals.

Older methods pool all BAMs when calling variants. This is necessary because a single low-coverage sample does not have enough data to recover hidden INDELs. However, this strategy is not that easy to massively parallelized; adding a new sample triggers re-calling, which is very expensive as well. As we are mostly doing high-coverage sequencing these days, the old problem with INDEL calling does not matter now. GATK has this new single-sample calling pipeline where you combine per-sample gVCFs at a later stage. Such sample combining strategy is perhaps the only sensible solution when you are dealing with 100k samples.

The so-called haplotype based variant calling is a separate question. This type of approach helps to call INDELs, but is not of much relevance to multi-sample calling. Also, of the three variant callers in your question, only GATK (and Scalpel which you have not mentioned) use assembly at large. Freebayes does not. Platypus does but only to a limited extent and does not work well in practice.

I guess what you really want to talk about is imputation based calling. This approach further improves sensitivity with LD. With enough samples, you can measure the LD between two positions. Suppose at position 1000, you see one REF read and no ALT reads; at position 1500, you see one REF read and two ALT reads. You would not call any SNPs at position 1000 even given multiple samples. However, when you know the two positions are strongly linked and the dominant haplotypes are REF-REF and ALT-ALT, you know the sample under investigation is likely to have a missing ALT allele. LD transfers signals across sites and enhances the power to make correct genotyping calls. Nonetheless, as we are mostly doing high-coverage sequencing nowadays, imputation based methods only have a minor effect and are rarely applied.


I think that genome-wide probabilities (i.e. not 0.25 for each base) should be fine. For a given segment of DNA, the chance that the sequence doesn't appear anywhere else in the genome is fairly high.
There is lots of patchy non-randomness (e.g. at centromeres, STRs, and other tandem repeat regions), but there's not much point in modelling that because it's not random.


Books on application of Machine Learning in Bioinformatics

One could read Bioinformatics: The machine learning approach as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...
...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.
I would like however to stress the importance of the basic background - one regularly sees announcements about breakthroughs achieved by reinventing techniques that were already described in the book that I mentioned.


Downloading a reference Genome for Bowtie2

It’s a matter of preference I guess but I recommend the Ensembl builds. Decide whether you want the toplevel or primary assembly, and whether you want soft-masked, repeat-masked or unmasked files. The naming schema is very straightforward; the combinations are described in the README file, and all files reside in one directory.

For example, if you want the unmasked primary assembly, the file to download would be Homo_sapiens.GRCh37.75.dna.primary_assembly.fa.gz.

As for GoldenPath/UCSC, there’s no need to download and concatenate separate chromosomes (contrary to what the other answer said); you can download the whole (toplevel) reference from the bigZips directory; from the README:

  This directory contains the Feb. 2009 assembly of the human genome (hg19,
  GRCh37 Genome Reference Consortium Human Reference 37 (GCA_000001405.1)),
  as well as repeat annotations and GenBank sequences.

There are essentially three options here:

chromFa.tar.gz, which contains the whole genome in one chromosome per file;
chromFaMasked.tar.gz, the same with repeats masked by N;
hg19.2bit, which is the whole genome in one file, but needs to be extracted using the utility program twoBitToFa, which needs to be downloaded separately.

In any case, I always download the reference and build my own index for mapping, since this allows me more control; not everybody might need this much control, but then building the index once is fairly fast anyway.


How to add OVERLAP_BP field to VEP output?

I also opened this as a Github issue, and received the following response:

  It seems there are some inconsistencies on our end. I'm looking into it.
  Until the fix, using --fields OverlapBP,OverlapPC works with any output format, while --fields OVERLAP_BP,OVERLAP_PC works only with default output format.

  A small note that this overlap data is available for structural variants only.


Software recommendations - DNA composition

There are more R packages available that calculate GC content, for example Ape's GC.content() function.

For example:

> library(ape)
> data(woodmouse)
> GC.content(woodmouse)
[1] 0.3873347

With a sliding window is in Biostrings package.

> library(Biostrings)
> DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC")
> letterFrequencyInSlidingView(DNA, view.width=4,letters="CG")

For codon usage bias, I think cubfits or sscu packages might be useful.


Will using smaller kmers help get larger contigs? If not, then what?

You can try s-aligner (free for 15 days) which usually gets quite larger contigs and quite larger NG50 for metagenomic data containing viruses (also phages). Indeed metaSPAdes also gets larger contigs, despite shorter than s-aligner.

Ref. s-aligner: a greedy algorithm for non-greedy de novo genome assembly


Appropriate value to replace outliers in clinical studies

In my experience, it is more common to use all clinical data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data.

If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.


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

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


You can use ifelse():

anno_maf$COSMIC <- ifelse(anno_maf$Hugo_Symbol %in% your_gene_list,
                          yes = "yes",
                          no = "no")


Why does the FASTA sequence for coronavirus look like DNA, not RNA?

That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil).

The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse transcribed, either by targeted reverse transcription or random amplification and thus converted to cDNA. cDNA is stable and is essentially reverse transcribed RNA.

The cDNA is either sequenced directly or further amplified by PCR and then sequenced. Hence the sequence we observe is the cDNA rather than RNA, thus we observe thymine rather than uracil and that is how it is reported.


How do I compare measured ratios to calculated ratios in peptide mass spectrometry?

From a statistical perspective ratios are icky because the variance of the sampling distribution of a ratio decreases with the size of the numerator and denominator being sampled to create the ratio. If you can get at the numerator and denominator separately, then you should absolutely do that and look at count based tests like the Chi-square test, Fisher's exact test, or Poisson regression.

That said, a t-test would be appropriate in the two-sample case if its assumptions are met.  You need to check that both samples are normally distributed, ideally via a qq-plot. You also need to check the variance of both samples and choose a variant of the t-test that does or does not assume equal variances.

If the samples are not normally distributed, then the next approach would be a Mann-Whitney U test.


Batch detection of CRISP proteins in fasta file

Is your goal to detect the actual proteins from the sequences? If you'll settle for detection of transcripts from the genes encoding these proteins as opposed to detection of the proteins themselves, here is an approach.

Aligning your transcriptome to a reference genome (probably with STAR) will be the first step. You have lots of short reads in your FASTA, but in their current state they are no more informative than a box of puzzle pieces. Aligning your reads to the reference will allow you to extrapolate which genes that you short reads correspond with. From here you can quantify the features (genes) in your transcriptome.

Lucky for you, there are only 3 genes in the cysteine rich secretory protein family, CRISP1, CRISP2, and CRISP3. Alignment with a reference (and annotations) and quantification of genes that your reads align to will give you all you need. You'll just have to find the three genes in your final counts table.

I recommend STAR as an aligner with the --quantMode GeneCounts argument. It's very quick if you've got the resources and is geared towards RNA-seq reads.


how to get access to genotypes and phenotypes used for a GWAS

How about:

dbGaP - https://www.ncbi.nlm.nih.gov/gap/ - in case you are searching for supporting datasets I leave its entry in BiŌkeanós

https://www.ddbj.nig.ac.jp/jga/ - the access might be a bit more complicated.

https://www.internationalgenome.org, not particularly large, in their FAQ they state that the individuals sequenced have not self-reported any diseases - perhaps you might want to add those with to your negative set.  Again, its entry in BiŌkeanós.

If you are looking for more comprehensive information, you could try biobanks - e.g. https://www.ukbiobank.ac.uk. These are usually open (or reasonably priced) for research but the approval process can take time.
A side note - some of GWAS catalogues will have a link to the study where the association originated from, e.g. here. You can sometimes find the datasource there.


Is there any tools that can generating FBA model from genome sequences or proteome sequences?

To generate flux balance, one should now the stoichiometric numbers of all the reactions that are taking place, as well as the conditions (that could affect the reactions). I am not sure we have enough knowledge, time and computationally power to do this.

At the moment I'm only aware of a paper simulating the Mycoplasma genitalium cell. But I'm not sure it was worth the effort as I haven't seen replicated or redone with other species, perhaps the model was not accurate enough to be worth it.

The tools used are here. The licenses are not a standard, but at least they seem free.


Exporting nt's differences from MEGA alignment

You don't say which version of MEGA you are using, but this works for MEGA X on Linux:

Open your alignment with the "Sequence data explorer"
At the "Highlight" menu, select "Variable sites".
Then at the menu "Statistics", select "Use only highlighted sites" and select what you consider an appropriate output - I choose "Display results as CSV".
Now at the menu "Data" select "Export data"; at the "Sequence data export options window" if necessary select the data format again, and select "For each site" at "Writing site numbers", and select "Only highlighted sites" at "Selected sites to include".

This will result in a comma-separated values file, showing only positions with variable sites.


Affymetrix tags (same ID's) present in different places of the genome

Is the annotation able to distinguish such 'same' tags, during DGE analysis?

No, you'll end up discarding such probes (assuming you have a reason to actually use microarrays still).

  Why aren't tags unique?

Platforms such as HG-U133 date back to 2002 (meaning they were designed a few years before that), when the human reference genome (and the reference transcriptome) were not nearly as "finished" as they are today. Probes that may have made sense given what we knew almost 20 years ago may no longer makes sense now given what we know. Further, microarrays are rather conservative, in that they tend to primarily gain probes over time and not have them removed (lest companies lose customers who are relying on given probes).


confidence ellipses for MDS plot in edgeR?

I can see that there would be a three step process to doing this:

Merge counts from all samples in the group and then resample pseudo-replicates from this. If x is a matrix with samples in the group being columns and genes being rows

s <- sample(row.names(x), n = mean(colSums(x)), probs=rowSums(x)/sum(rowSums(x))
stab <- table(s)
s <- as.vector(stab)
names(s) <- names(stab)

Do this thousands of times.
You would then want to project those onto your MDS plot space - i'm less clear about how you could do this without perturbing the space itself.
Calculate the parameter for an ellipise that would contain 95% of these points. Again, I'm not so sure about how to do this, but I think the car pacakge might be a good place to start looking.


What is the state of the art for GWAS in terms of the statistical algorithm for either Case/control and Quantitative traits?

I've got the answer to my question at https://www.biostars.org/p/410466/ by chrchang522

"The main regression executed by Plink was introduced by EIGENSTRAT in ~2006; see https://www.nature.com/articles/ng1847 . This is actually straightforward to write in R/Python from scratch; the harder part is optimizing the implementation for large datasets.
The Firth regression added to Plink 2.0 to improve handling of rare variants and imbalanced binary phenotypes was motivated by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4049324/ .
Mixed linear models provide better statistical power when you have lots of close relatives in your dataset, but are much trickier to solve; actually, this is still a significant research area. Two tools covering parts of the current state-of-the-art are SAIGE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6119127/ ; handles imbalanced binary phenotypes, but relatively slow) and fastGWA (https://www.nature.com/articles/s41588-019-0530-8 ; great speed, but doesn't support dosage data yet and uses a misspecified model for binary phenotypes)."


How to use busco installed through conda?

From the README, it looks like the command-line entry point is not busco but run_BUSCO.py. This seems to agree with the bioconda recipe, where run_BUSCO.py and run_busco are listed among the commands.


How many single copy gene are needed to construct the phylogenetic tree?

The one thing that strikes me is that OrthoMCL was written under David Roos supervision, which really means it was written with Toxoplasma gondii in mind, i.e. a eukaryotic parasite with notable levels of gene duplication. They tested it against Plasmodium falciparum, but that isn't surprising because both are members of the Apicomplexa. Personally I would suggest looking at a tool more orientated towards bacteria, albeit I'd like to see the output of OrthoMCL.

Prokka is a tool that might give better results, here, which is designed for bacteria, but personally I would extract the orthologes/homologues directly if the primary purpose is to contruct a phylogeny, i.e. I would parse annotated genes which are not "hypothetical".


Least present k-mers in the human genome

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

From the instructions in the user guide:

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

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

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

jellyfish histo mer_counts.jf

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

jellyfish query mer_counts.jf AACGTTG

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

jellyfish dump mer_counts.jf > mer_counts_dumps.fa


Large-scale gRNA design for a CRISPR screen

Perhaps GenomeCRISPR, I've personally never used it, but it has a well-documented API and appears you could automate it to go through the whole genome.


Whole genome aliment can be done using Progressive Mauve, LAST or Mummer. For bacteria I used Mauve since it has also very nice visualisation engine. A very new tool is Minimap2, a super fast mapper that supposed beside read mapping be able to handle reference vs reference. However, I do not know how performance of it compares to the tools mentioned above.

If you are interested in a rough idea of the shared genome regions, you can use bevel. Bevel is not really an aligner, it is more like a dot-plot, but it is super fast (even for mammalian sized genomes).


Is there a database of protein sequences/structure along with their melting temperature?

ProTherm2 contains thermodynamic properties of proteins and mutants. https://www.iitm.ac.in/bioinfo/ProTherm2/ .
NAR Database issue advertises it as containing melting temperatures.
For the reference - found it with https://biokeanos.com/search?q=protein+%2Bmelting


Mapping statistics from bam file using bbtools and sambamba

Try samtools flagstat or sambamba flagstat to get similar type of information about your bam file. It may be easier to interpret. Both of these tools can be installed using conda if needed.
EDIT:
Yes, the results of sambamba flagstat you showed in your edited post indicate exactly what you said: 94.91% mapped and 100% - 94.91% = 5.09% unmapped:
64035749 + 0 mapped (94.91%:N/A)


Customizing colors in geom_arc_bar (ggforce)

For a discrete scale your Code 1 example was nearly there.  Just use scale_fill_manual instead of scale_color_manual.  For more control over how the discrete values map to color values you can name the vector, too.

my_colors <- c("red","white","blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)
# or, if you want consistent output independent of the order of
# df$discrete_fill:
my_colors <- c("1" = "red", "2" = "white", "dummy" = "blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)

You can also just use a single call to geom_arc_bar for the whole data frame (rather than a list of call results from each row) and you'll get the same effect.

ggplot(df) +
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie") +
    theme_no_axes() + scale_fill_manual(values=my_colors)


multiple COSMIC id for the same mutation

Harry Boutselakis addressed the question via the customer support of the COSMIC web site.

  Dear Ivan, Thank you for using COSMIC.

    I would like to know why there are different entries for the same mutation   in COSMIC, for example COSM6954941 and COSM12833 both refer
    to ERBB4   c.908C>A.

    In this specific case the field Gene name in the two entries is different,   COSM6954941 report ERBB4_ENST00000402597 while COSM12833
    report ERBB4,

  The mutations are annotated onto two different isoforms of ERBB4. So
  this annotation refers to two different transcripts. In these cases
  the consequence of the mutation on the mRNA and polypeptide is the
  same reflected by the syntax.

    but there are many cases in which the Gene name is the same, for example   COSM1361 and COSM1219.

  This case the consequence one mutation on the mRNA is unknown , but
  the amino acid change is the same on both cases.

  The information stored in COSMIC is curated from experimental data,
  depending on the evidence and experimental methodology used, the
  mutations are differentiated accordingly. I hope this answers your
  question.

  Kind regards,

  Harry


Mapping a list of cells in seurat featureplot

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

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

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

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

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


What is the difference between SAM mapping quality and Blast E-value?

The E-value and the mapping qualities are two very different things.

The E-value is "a parameter that describes the number of hits one can 'expect' to see by chance when searching a database of a particular size". More details can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect

The mapping quality is an attempt to estimate the probability that a given base from a given read is mapped correctly to a particular place in the reference genome. Different aligners calculate mapping qualities very differently, so there is no simple way to compare them across aligners. But in general the mapping quality will consider things like: (i) the quality of the base call; (ii) the repeat structure of the reference; (iii) the alignment algorithm; (iv) whether the read has a mapped pair; (v) anything else the author of the mapping software thought might help.

E-values and mapping qualities are, therefore, measuring two fundamentally different properties with very different uses.


Mapping bisulfite reads to a reduced size genome

Don't restrict the portion of the genome you're mapping to, you'll just increase false-positive alignments by doing so. If you're worried about mapping time, then (A) use more cores or (B) use a different mapper (BSMAP was the quickest out there last I looked, though in my opinion it's fine to wait a bit longer for better quality results (I normally suggest bwameth to people)).

As pointed out by user 719016 in the comments, the above does not apply to amplicon sequencing. With that, you'll get VASTLY better performance if align to a reduced-size genome...and without an appreciably increased false-positive rate (unless your target-enrichment for sequencing is poor).


Calculating L2 norm (Euclidean norm) to set of genes in time series

If you read your matrix into R and saved it as m:

diffMatrix = m[,2:9] - m[,1:8]  # The change per-timepoint
l2n = apply(diffMatrix, 2, function(x) sqrt(sum(x^2)))  # The L2-norm

R also has a norm() function, which can directly produce the L2-norm, but it turns out that this is slower than the simple method shown here. If you click through to the other question you can find somewhat faster implementations of this, but anything will suffice for small datasets like yours.


Defining methylation level of CpG island

In my experience, the most common approach to find DMR is to use the DSS R package to call and merge DML into regions.
DMC/DML (Differential Methylation Cytosine/Loci) is often calculated by comparing the counts of methylated and unmethylated cytosine at one position between samples.
Steps are to first call DMLs with the callDML() function. Then call DMRs from these first results with callDMR(). Finally annotate the DMLs or DMRs using a GFF file. From there you can look for DMR falling into CpG island. Some of the parameters that can be tuned are the minimum region length, minimum number of differential CpG sites, percentage of CpG sites required to be significant. However I have never heard of what would be the standard values for these parameters. It may not be clearly established yet.
So because the process of calling DMLs and DMRs need robust statistical testing, it is suggested to first call the differential methylated states and then annotate it.


Uniprot: Determining protein and organism given a DNA sequence

Well, let's see. The first step is to simply run a nucleotide blast (blastn) against the nr (non-redundant) database. I did this using the NCBI blast tool. While, as expected, I got many hits, there was only one hit which matched the entire length of the query sequence perfectly:

Diplonema sp. ATCC 50224 cytochrome c oxidase subunit 1 (cox1) mRNA, complete cds; mitochondrial
Sequence ID: JF698657.1
Length: 1612Number of Matches: 1
Range 1: 141 to 770
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
1164 bits(630)  0.0 630/630(100%)   0/630(0%)   Plus/Plus
Query  1    ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  60
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  141  ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  200

Query  61   GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  120
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  201  GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  260

Query  121  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  180
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  261  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  320

Query  181  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  240
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  321  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  380

Query  241  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  300
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  381  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  440

Query  301  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  360
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  441  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  500

Query  361  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  420
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  501  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  560

Query  421  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  480
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  561  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  620

Query  481  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  540
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  621  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  680

Query  541  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  600
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  681  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  740

Query  601  CATCATACTTCCAGTATTCGGTATAGTCAG  630
            ||||||||||||||||||||||||||||||
Sbjct  741  CATCATACTTCCAGTATTCGGTATAGTCAG  770

So, we now know the sequence is from the gene for mitochondrial Cytochrome C Oxidase subunit 1 from Diplonema sp. ATCC 50224. That the sequence is mitochondrial, is very important since the mitochondrion of this species has a slightly different genetic code. If you go to the NCBI's page for this sequence (JF698657.1), you will see this line:

/transl_table=4

That tells you it is using translation table 4:

4. The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code (transl_table=4)

    AAs  = FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
  Starts = --MM------**-------M------------MMMM---------------M------------
  Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG

See how all of your translated reading frames have in-frame stop codons? That's because you translated them using the standard genetic code instead of the correct mitochondrial one. As you can see above, in this code, the codon TGA is not a STOP, as in the standard code, and instead codes for tyrosine (W).

So, if you look at the result of the protein blast you ran on UniProt, you will see:

Query1   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL   59
         LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI FL
F5CI61   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL   107

Query60  VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   118
         VCGLLLITGSSWMEEGSGTA TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
F5CI61   VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   167

Query119 CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   177
         CTIAYARRSLMAMLQSSLYP SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
F5CI61   CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   227

Query178 GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  206
         GDPMLYQHLFWVFGHPEVYIIILPVFGIV
F5CI61   GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  256

The only locations where there isn't a 100% identical match is where a W is aligned to a STOP (- in the Uniprot output). So, your actual protein sequence is the 2nd frame of the 5' to 3' direction which, when translated using the correct code, is:

LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL
VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV


RSeQC provided bedfile does not work with read_distribution.py

Ok, I figured it out already, it's and old issue I come accross every now and then. The original bed file lists the chromosomes as chr1,chr2,...,chrX whereas my bam files uses 1,2,...,X

I fixed the bed file like this (python3):

import pandas as pd
bedFile = '/hg38_Gencode_V28.bed'
bedFile = pd.read_csv(bedFile, sep='\t', header=None)
bedFile[0] = [value.strip('chr') for value in bedFile[0]]
newBedFile = 'hg38_Gencode_V28_reformatted.bed'
bedFile.to_csv(newBedFile, sep='\t', header=None, index=None)

And now read_distribution.py works.


Inspection of gene expression in scRNA-seq data

In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the point the course instructors are trying to get across with that message box.

Eyeballing the results, you should also exclude any ERCC genes as well, as they are synthetic constructs and not from the target samples (although record their counts, as it may be informative for QC/normalisation purposes); they'll have odd distributions because of this.

ERCC sequences are typically added based on the concentration of DNA in a sample, rather than the number of cells. However with single-cell data (i.e. with known cell counts), it'd be possible to do a proper ERCC-based normalisation, assuming an equal amount were added for each cell. I suspect that knowing a bit more about the ERCC addition protocol would help in understanding why their distribution is so different.


How to map short sequences to long reads, recovering all multiply-mapped high-quality matches

LAST has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see here), as well as abundant base calling differences present in single-molecule sequencing.

Minimising error is not necessarily going to be the best option, and concentrating only on unique hits will miss a lot of signal. There are frequently multiple identically plausible positions, even at zero error, and the long reads will have their own associated errors.

I also find the limitations a little odd. Canu can do read-level correction on long reads based on overlapping reads, and if assembly is possible, then nanopolish can be used to correct nanopore reads for systematic base-calling error introduced by the methylation of unamplified template.


How are different forcefields developed/evaluated?

Some force-fields rely solely on theoretical physical principles, while others are tweaked specifically for the goal of accurately depicting certain structures like sheets and helices or movement/diffusion aspects at different theoretical temperatures of a set of well-studied proteins. (Top-Down vs. Bottom-Up)
I haven't done MD myself but I am in an MD group and people seem to normally stick to a single force-field for a given MD project but it happens that they do very short simulation runs with force-fields just as a test.
They definitely don't plan on doing the full simulation on different force-fields unless the comparison of those is their main goal.


Minimap2 -ax map-pb doesn't output tlen field

The TLEN for single-end reads is defined as 0 in the SAM specification, which is why this happens. You can still get the end position of the alignment by parsing the CIGAR string. An example using the pysam package in python is below:

import pysam

bam = pysam.AlignmentFile("myfile.bam")
for b in bam.fetch():
    print("This alignment spans {}:{}-{}".format(b.reference_name, b.reference_start+1, b.reference_end))

The reference_end value is automatically calculated from the CIGAR string.


why in RNA seq don't we only use reference transcriptome?

Pseudoaligmnent to a reference transcriptome with something like Kallisto is perfectly valid, assuming you are in an organism whose reference transcriptome is very well documented.
But if you want to find novel transcripts or novel splicing junctions, you need to align to genome.


Pre-filtering genes for Principal Component Analysis

You can filter out genes that do not show variation across samples, they would not differentiate samples anyway. For the specifics on how to do so, please see the Subpopulation Analysis section of a nice single cell RNA-seq workflow from the Kharchenko Lab. I believe there won't be fundamental differences when applying this approach to bulk RNA-seq.


Run Rscript for multiple samples (files)

Please do not use setwd in a script, it changes working folder without explicit specifications from the user, so the results could be unpredictable for people that are not familiar with the content of the script.

Create a tab-separated file with input filename and output filename, or at least a file with all input files, unless you can pick them using a list.files(pattern = "<pattern_here>"). Read this file as a data.frame and use it with an apply, so a function processes the df row-by-row. Let's say the file (named metadata.txt) is:

in_file out_file
tgf.txt tgf.output.txt
abc.txt abc.filtered.txt
xyz.txt xyz.randomname.txt

Define a function called process_row() in the Rscript:

process_row <- function(r) {
    d <- read.table(r[[1]], header = FALSE);
    x <- aggregate(d[, 2], list(d$V3), median);
    write.table(x, r[[2]], quote=FALSE, row.names=TRUE);
}

And do the actual work:

meta_data <- read.table('./metadata.txt', sep = ' ', header = TRUE, stringsAsFactors = FALSE);
apply(meta_data, 1, process_row(r))
# You don't need a q() in an Rscript because it exits after the last line of code anyway.


I'll stick my neck out and simply state that allopatry alone would account for the pattern you are observing. Thus populations can exihibit variation in allelic profiles, particularly where Mendellian inheritance is operating (i.e. higher eukaryotes), singly due to geographic isolation.
This represents the Meyers view of population genetics and would account for variation within neutral theory, e.g. silent mutations. Obviously I'm guessing the populations were selected by geographic foci in generating the allelic frequencies.


What do negative dN or dS mean in cal_dn_ds in Biopython?

dS = -1
.. there are two suggestions

I suggest its an error, possibly based on nucleotide calls within the S, check for anomylous nucleotides, particularly over a non-synonymous codon
It might mean their are no synonymous mutations. It is not the best idea to assign this a value of 0 in a ratio because it will cause the ratio to fall over when in fact their could be strong selection pressure, ie non-synon mutations. Having an incrementally small value is problematic because dN if positive could anything up to infinity, so -1 might flag this problem.

Okay I now know why... You are correct to be concerned with the result because its serious.

My sincere advice to you is do not use any dN/dS ratio where -1 is occurring under any circumstances. On this occasion I'm not going to explain why for two reasons:

First and foremost, as a phylogeneticist you should know why, it relates to the fundamentals of the discipline.
(Hint) The answer is in my past posts on this forum.

To conclude, I'm pretty nuanced in statistical dogma, you can check my record, IMO most dogma is tribalism because every calculation has weaknesses. HERE HOWEVER, there is a real problem and again I urge you to ditch any data set where -1 is occurring and everything associated with it. That is a fairly dire warning should give you a clue about whats wrong. If you are not prepared to do this, then you must report the result and the situations where it is arising (don't hide it).


samtools / bamUtil | Meaning of  as Reference Name

As per the SAM/BAM format specifications:

  "RNAME: Reference sequence NAME of the alignment. If @SQheader lines are present, RNAME(if not ‘*’) must be present in one of the SQ-SN tag. An unmapped segment without coordinate has a ‘*’ at this field.  However, an unmapped segment may also have an ordinary coordinate such that it can be placed at a desired position after sorting.  If RNAME is ‘*’, no assumptions can be made about POS and CIGAR." -- http://samtools.github.io/hts-specs/SAMv1.pdf, emphasis mine

So, yes, these are the unmapped reads.

I haven't used this splitChromosome utility before myself, but as it splits the BAM file into one BAM file per chromosome, I suspect it simply also creates a BAM file for the unknown (*) chromosome. Why the numbers don't match I can't tell.

How to deal with that file in your pipeline will depend on what the pipeline is meant to do. But if you aren't interested in these reads, I suppose you could simply remove it before continuing to the next step.


How to combine multiple kernels of large sample datasets?

Welcome to Bioinformatics Stackexchange @SD1024. A recently published algorithm known as MOFA (Multi-Omics Factor Analysis, paper, github) is generating a lot of interest, and is designed answer exactly the sort of biological question you are describing. It claims to extract axes from multiple matrices with overlapping samples, but not necessarily overlapping features and handle missing data elegantly. I havn't tried it myself you, but I'm quite excited to do so.


I'm not sure Seurat is the best tool for this as it was developed for single cell RNA seq data and there are a few intricacies of that type of data that are very different from bulk RNA seq.
For bulk there are really good packages available and corresponding workflows, e.g. limma, edgeR and DESeq2.

The main problems with Seurat for bulk RNA-seq:

Seurat expects counts as input - FPKM are not counts nor are they log counts or log norm counts. It's a specific normalisation method that takes into account gene length and library size and breaks the link between gene counts and variance.
Because of the distribution of reads from a single cell that you get from the sequencer Seurat scales all unique reads (i.e. UMIs) to 10,000 per cell and usually regresses out the nUMI count. Again, becuase you've got FPKMs this makes no sense.
The other functions you can use on bulk data, i.e. tSNE, finding clusters and markers (which is essentially DGE analysis). The one thing to note for the latter is that again because you've got FPKMs they don't follow the negative binomial distribution. Your best bet is limma trend but I'm not sure it's included in Seurat.


How do I export my biom object to a file in biom format?

You can use the functions to_hdf5 and to_json from the biom-format package:

to_hdf5

with biom.util.biom_open('output.biom', 'w') as f:
    merged_table.to_hdf5(f, "example")

to_json

with open('output.txt', 'w') as f:
    f.write(merged_table.to_json("example"))


What are the coordinates of the extended major histocompatibility complex (xMHC) in GRCh38?

Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13
which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.
The "extended MHC" defined here [1] is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The xMHC is defined in [1] as the regions between HIST1H2AA and RPL12P1 genes, which would give you chr6:25726063-33400644 in Grch38. This is very similar to your HIST1H2AA - KIFC1 definition (chr6:25726063-33410226). Of course if you're trying to reproduce a result, probably better to use the original paper's definition.
[1] https://www.nature.com/articles/nrg1489


Gene function annotation - bacterial genome

The PATRICBrc.org is a free resource you can annotate your genome and compare it with publicly available genomes.

Refer this tutorial for more detail.

Once you have annotated your genome in PATRIC, you will have PATRIC Local families (k-mer based protein family), which will help you identify all genes with a specific function.


Solutions for managing data in a small bioinformatics / 'omics lab?

I don't know of any prebuilt products, but I can describe how we managed this in my postdoc lab, and how I plan to manage it in my newly-started group.

Rule 1: All work happens in the projects project directory on the central
filestore, not on your desktop or laptop.

Rule 2: Heavy computational work is done by the groups standard analysis pipelines. Interpretation is done in jupyter notebooks or Rmarkdown.

A project has a directory on the group's filestore. That directory has a fixed structure:

proj001----raw_data             *
       |
        --external_datasets
       |
        --src                  *
       |
        --notebooks            *
       |
        --pipeline1
       |
        --pipeline2
       |
        -- etc
 * link to a seperate, backed up filesystem.

Pipelines are where heavy duty analysis happens and everyone uses the same standard pipelines, producing the same standard output files and database structure.

So common pipelines might be the mapping pipeline, the readqc pipeline, the differential expression pipeline, the exome pipeline, the chip-seq pipeline etc...

Pipelines have three important outputs: an automatically generated html report, an sqlite database and files in the export directory.

When we archive a project, this is what we save, along with the pipeline's configuration file and log file.

So if I know that Jane Bloggs did an RNAseq 5 years ago on a cell type i'm interested in, if I konw that that was project 5, I know that in the project 5 directory there will be a diff_expression_pipeline directory, and that it will contain an sqlite database called csvdb and that that will have a table called refcoding_deseq_gene_diff and that table will follow a known format. There will be bigwigs in the export directory. Or the BAMs will be in the mapping_pipeline directory.

Of course the problem remains knowing that Jane Bloggs did this RNAseq and what the project was called. We use a Wiki for this, but its not ideal.


Is there a standard method to validate whether a bigWig file is correct?

There's no standard validation method, but a simple one would be python -c "import pyBigWig; bw = pyBigWig.open('your_file.bw'); print(bw.IsBigWig());". If you get True then it's a bigWig file. If you get False or an error then it's not (False indicates a bigBed file).

Alternatively, file your_file.bw should return data.


Does BLASTClust guarantee that proteins in different clusters are dissimilar?

Based on the documentation (that I linked to in comment), it seems that the clustering algorithm used by blastclust is single-linkage hierarchical clustering.
SO what you can expect is that A will be placed with the protein to which it is most similar. And then the last protein will be added as an outgroup to that cluster due to its similarity to A (assuming that there are no other proteins).
So if your goal is "I want clusters that are guaranteed to contain dissimilar proteins", then you might not be able to naively take blastclust output and interpret it as such. You could either postprocess it to find the set of clusters in which no protein is more than 30% similar to a protein in any other cluster (which may just be one big cluster!), or find another clustering algorithm that is more suited to your task.
Update
It does look like blastclust uses a minimum similarity (not identity) threshold (from blastclust docs). So they use an implementation of single-linkage clustering that stops before merging everything into a single cluster:
If the coverage is above a certain threshold
 AND
the score density is above a certain threshold,

these two sequences are considered to be neighbored.

Thus determined neighbor relationships is considered symmetric and provides
the base for clustering by a single-linkage method (which puts a sequence
to a cluster if the sequence is a neighbor to at least one sequence in the
cluster)."

I don't think that this is a general property of the generic clustering algorithm though (see the "Naive algorithm" section).


How to run enrichment analysis of protein functional annotation?

I would recommend R instead of sticking to websites.

There are many tools for enrichment analysis. I like to use GOseq, which is initially made for RNAseq data, but can also be used for proteins (I have used it for proteomics). You can use the length of the protein instead of transcript length (to correct for length bias), or you can exclude the length information by just using the hypergeometric method.

See manual for code.

Subsequently you can use gogadget, a tool I made for visualization of GOseq results with a heatmap or cytoscape network.

Additionally, you might wanna take a look at clusterProfiler. With this package you can do enrichment analysis for many clusters in parallel.


How do you normalise read coverage in a BAM file?

I just updated one of my SAM-normalisation tools to more accurately normalise mapped reads to a target coverage. It uses reservoir sampling on a SAM-formatted file to choose reads to bring the coverage up to 100X (if possible).

https://gitlab.com/gringer/bioinfscripts/blob/master/samNormalise.pl

Example Usage:

samtools view -h input.bam | ~/scripts/samNormalise.pl | samtools view -b > normalised.bam

The output BAM file can be converted into a fastq file using 'samtools fastq':

samtools fastq normalised.bam > normalised.fastq

Here are some example before and after shots ['samtools mpileup input.bam'; 'samtools mpileup normalised.bam']:

Notice that where the coverage is below 100, the coverage is identical. Coverage is below 100 in some areas because some reads end, and no additional reads start until later on in the reference sequence.

The options are coded as variables at the beginning of the script at the moment, but I can easily convert them to command-line options in the future if there's demand for it.


Error while using E-utilities on slurm HPC

501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed) indicates that a specific perl module is missing.
If indeed your edirect installation is from March 2019 (as I interpret from e-utilities/03.02.19), you should update your edirect installation to the latest version. If you update, I expect the necessary perl module will be installed and the error will go away.


how to calculate coverage of each base of mapped reads?

You are looking for bedtools genomecov, specifically -d option. See https://bedtools.readthedocs.io/en/latest/content/tools/genomecov.html , for example:

bedtools genomecov -i A.bed -g my.genome -d

I recommend not reinventing the wheel: use bedtools instead of reimplementing the same functionality in Perl.


This question will be a little hard to answer without more information.
For example, we will need to see your dataset (whole_root_pred), to decide why Stunting_Root is NULL.

You might need to initialize Stunting_Root as a variable. It is currently not clear if it is e.g. a column of your dataframe, or just uninitialized. Uninitialized variables are NULL, which would explain your problem. randomForest might not know to look for strata inside your dataframe, for example. Is it in your dataframe?
Also, I might be missing something, but why are you passing ? as a response? I'm not an expert but I believe that is an illegal character in R (I'm pretty sure?).


Why are there missing calls in a VCF file from exome sequencing?

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd filter out the VCF to only include on-target loci before proceeding with further analyses.

If you did your variant calling on samples separately and then merged them together, you'll have missing calls at loci called variant in one sample and homozygous reference in the other.

You can get over this by asking the variant caller to emit all sites within the target region, which would lead to a large VCF or do multi-sample calling, which is more computationally intensive.

One other way is to do variant calling on all samples separately, generating a BED file of all variant regions in the cohort and then asking the variant caller to genotype those loci across all samples before merging them together.


Hierarchical clustering for outlier detection - single cell RNASeq & WGCNA

For outlier identification I suggest using the sample network approach developed by Oldham et al. It basically amounts to constructing a inter-sample connectivity (assuming normalizedData contains the log-normalized data with genes in colunms and samples in rows)

k = colSums(cor(t(normalizedData))),

scaling the connectivities

Z.k = scale(k)

and then removing samples with high negative Z.k (e.g. Z.k < -4).

A slightly different approach is based on euclidean distances rather than correlations,

Z.k = scale(as.matrix(dist(normalizedData)))

and removing samples with high positive Z.k (since distance is a measure of dissimilarity).

Regarding the scale-free topology fit plot, I wouldn't worry about it too much, the fit index sometimes jumps around a bit. It certainly does not look like you have strong outliers or very heterogeneous data. But it is always a good idea to plot the sample clustering tree to check (1) that there are no strong clusters, and (2) that the outlier identification did something reasonable.


How BLOSUM Matrix is constructed and calculated

Sean Eddy wrote an excellent explanation of the BLOSUM Matrix in his Nature Biotechnology paper Where did the BLOSUM62 alignment score matrix come from? He also gives real examples of how one would calculate the score for different amino acid pairs and where the numbers come from.  Further references given at the end of the paper might also be useful reading, specifically: Amino acid substitution matrices from an information theoretic perspective, Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes, and Amino acid substitution matrices from protein blocks.


A program "PSICOV" "stopped" from a few days on Ubuntu 18.04?

In summary just be patient. psicov doesn't  parrallelise so a tenfold increase in CPU will not hurry it along.

Also psicov is not RAM intensive given the machine you are describing, so you haven't hit a RAM bottleneck.

I don't know the ps output format you are using I'm use to seeing %CPU and %MEM from ps -aux on Susie and OSX. You should also try,

top -c -u username

or for OSX top -u

You should see the psicov happily chuntering away consuming ~100% of a single process. There is a UNIX/Linux guru specialist on BioSE who could almost certainly provide better diagnostics.


That's a warning. Presumably, a match failed ($& is the default Perl variable which holds the result of a match operation). So something in your input is not exactly what the script was expecting. You can probably safely ignore these messages.

Alternatively, update your question and add a few example lines of input which we can use to reproduce the error and maybe we can figure something out. If that isn't possible, at least show us line 1962 of the annotate_variation.pl script since that's the one that's failing.


How to convert BAM file to bigWig in Python?

deepTools has a (somewhat poorly documented) API, since it's a python package too. The basic code framework is:

from deeptools import writeBedGraph
from deeptools.getScaleFactor import get_scale_factor

wr = writeBedGraph.WriteBedGraph(options and input files)
func_args = {'scaleFactor': get_scale_factor(...options...)}
wr.run(writeBedGraph.scaleCoverage, func_args, output_file_name, ...)


Splitting characters of a column in repeated lines and reducing them to one line

OK. This is a bit convoluted, but it works. There might be simpler solutions...

First we'll write a function that takes a data frame and a character position (ie; first or second position in the string), groups the data frame by the second and third columns, counts occurrences of each character in the specified position of the 18th column, and returns a vector according to your ranking scheme (this doesn't include counts of N since it doesn't occur in your example data, but it's a simple extension of the case_when line to add it):

library(dplyr)
library(tidyr)
library(stringr)

summarise_characters <- function(df, pos) {
  df %>% mutate(character=str_sub(X18, pos, pos)) %>%
  group_by(X2, X3, character) %>%
  summarise(n=n()) %>%
  ungroup() %>%
  spread(character, n, fill=0) %>%
  mutate(pos1=case_when(`@` > ß ~ '@',
                        `@` == ß ~ 'N',
                        `@` < ß ~ 'ß')) %>%
  `$`(pos1)
}

Next we remove duplicate lines from the data frame, then create a new column, pasting together the output of our function called on the first and second character:

df %>% group_by(X2, X3) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(counts=paste0(summarise_characters(df, 1), summarise_characters(df, 2))) %>%
  select(X1, X2, X3, X4, X11, X12, X13, counts)

# A tibble: 3 x 8
     X1        X2        X3                   X4       X11                  X12       X13 counts
  <int>     <int>     <int>                <chr>     <int>                <chr>     <int>  <chr>
1     1  90042011  90042031 AGTGCCACCAGGTGGGGCCG  90184354 GAGGCCACAAGAGGGCACAA  90184337     @N
2     1  94853396  94853416 GCCGCCACTAGATGGTGCTA  94969274 CAGTCCTCTAGAGGGAGCCC  94969266     ß@
3     2 103584283 103584303 TGGGCCACTAGGAGGCACTG 103841456 CATGCCACAAGAGGGCATCA 103841456     NN


Show presence of known mutation in RNA-seq data