What methods exist to calculate RNA expression profile similarity

There's a new paper that's just appeared on the subject of "proportionality", including a method by which RNA expression might be compared.

This is a new concept for me, and the article is not easy enough for me to read that I can write up a quick summary; the authors don't seem to devote a section in the paper to defining "proportionality". However, here's an interesting chunk from the article:

  We graphed the network of relationships between these mRNAs (S5 Fig.), an approach similar to gene co-expression network [12] or weighted gene co-expression analysis [13] but founded on proportionality and therefore valid for relative data. The network revealed one cluster of 96, and many other smaller clusters of mRNAs behaving proportionally across conditions.
  ...
  We are also keen to raise awareness that correlation (and other statistical methods that assume measurements come from real coordinate space) should not be applied to relative abundances. This is highly relevant to gene coexpression networks [12]. Correlation is at the heart of methods like Weighted Gene Co-expression Network Analysis [13] and heatmap visualization [14]. These methods are potentially misleading if applied to relative data.


We can plot using lists as well, using fromList function, see example:

library(UpSetR)

myGeneSets <- list(
  set1 = c("gene1","gene2","gene3"),
  set2 = c("gene1","gene4","gene5", "gene7"),
  set3 = c("gene1","gene7")
  )

upset(fromList(myGeneSets))

  fromList

  A function to convert a list of named vectors to a data frame
  compatible with UpSetR.

fromList(myGeneSets)
#   set1 set2 set3
# 1    1    1    1
# 2    1    0    0
# 3    1    0    0
# 4    0    1    0
# 5    0    1    0
# 6    0    1    1


Error in .checkedCall : subset indices out of range

Just posted the question and found out why the error was happening. I had two datasets: raw, and filtered one with half of the rows as raw. I did normalization on the filtered one. Then I calculated is.mito using the unfiltered, raw dataset. So, it produced indexes that were outside of the range of the dataset.


mrbayes, Sumt ; Error when setting parameter "Filename" (1)

Having thought this through, I think you've made a mistake. The file it wants is the outfile, you have given it the infile only. My original answer is below.

====================

From memory (my personal memory, not RAM, SSD or Stata), in the run command, you need to specify,

nruns=2

sumt will then combine both runs. If the file name is mystuff.run and nruns=2, there will be an output of mystuff.run1.t and mystuff.run2.t.

If nruns isn't set then sumt will not work, or if the MCMCMC output ends up with different names again sumt will not work.

With the burnin, you simply specify this in MrBayes. It is easy to write a script that will perform a 25% burnin. Tracer is nicely written, presumably C, its fast, but it has been around for a long time.


Downloading an open-source MAF with tumor and normal BAMs?

There's the NCI Genomic Data Commons:

  For all submitted sequence data, including BAM alignment files, the GDC generates new alignments in BAM format using the latest human reference genome GRCh38 with standard alignment pipelines. Using these standard alignments, the GDC generates high level derived data, including normal and tumor variant and mutation calls in VCF and MAF formats, and gene and miRNA expression and splice junction quantification data in TSV formats.

But that needs an NIH account for individual-level data, which might be too restrictive for you.

There's also the Texas Cancer Research Biobank (7 human cancer patients with matched normal):

  These tumor and normal specimen sequence files (FASTQ and BAM), somatic variant calls (VCF and MAF) and germline MAF files are available for each consented patient as described on the TCRB Open Access Privacy page.

That also requires account registration, but the registration process is automated and much simpler (requires just a name, and acceptance of usage conditions). There's a published paper about that dataset.

I had some vague feeling that another repository existed (created within the last year), but I can't seem to find it from quick web searches.


Error: package or namespace load failed for 'ChIPseeker'

Simple, use the latest version of R and Bioconductor, now on version 3.6.1 and Bioconductor 3.9. Those are the versions that are supported. If you can't upgrade R, check that BiocManager::valid() is true, you might have ended in some inconsistency between Bioconductor versions.

Looking at your sessionInfo, you seem to have some problems with the installed libraries, check and correct those too if possible


Protein folding problem model for quantum computers?

Protein folding problem can be viewed as a minimization problem. One can use quantum annealing to perform the minimization. Running this on quantum computers would improve the performance since they can perform the tunneling directly.

In fact, quantum annealing was used (blog post) for lattice protein folding on the the D-Wave quantum computer (128 qubits). But this was an extremely simplified model, with only 40 discrete states.


DMPfold: is it suitable to get 5000 models to compare with rosetta

DMPfold tends to produce low-variability models - it produces a small number of distinct conformations, and either the right one is in there or it isn't. This arises from the way the models are generated from distance constraints. Rosetta, by contrast, produces many different conformations which is why 5000 or more models tend to be generated. This arises from Rosetta's method of assembling small fragments, which allows for many more combinations.

You could generate 5000 models with DMPfold. Running with 10 5000 would give you 5000 models as output from the last iteration. Running with 10 500 would give you 500 models from the last iteration but you could also use the 500 models from each of the previous iterations (saved as ensemble.1.pdb etc.) to get 5000 models overall.

However I wouldn't recommend it, it sounds like a waste of disk space and compute time. Just run with default parameters 3 50 and use the final models saved as final_*.pdb.

For some numbers on how DMPfold compares to larger numbers of Rosetta models, and how much better DMPfold gets when you use more models, see Table 1 of the paper.


PacBio long-reads impact in transcriptome de novo assembly?

A few comments:
Never use N50 as a metric especially for transcriptomes. It has some semblance of relevance for genome assembly, but all that is void for a transcriptome with inherently dynamic lengths.
At the end of your IsoSeq pipeline, you should have (ideally) full-length transcripts. Have you considered that you don't miss much in your IsoSeq data? The biggest gain using IsoSeq is that it should have fewer non-existent isoforms assembled than a De-Bruin graph approach.
Let's say you have the following locus in your genome:
5'                                                                 3'
------------=======-----------=======----------=======-------------
Upstream      E1     Intron     E2    Intron     E3      Downstream

There are theoretically 7 isoforms:

E1
E2
E3
E1,E2
E2,E3
E1,E3
E1,E2,E3

A short read assembler has to look at evidence that connects exons to see which isoforms are represented. Let's say you have the following PE reads:
5'                                                                 3'
------------=======-----------=======----------=======-------------
               -->..............<--- Read 1
                                 -->...............<-- Read 2
                 --___________-->..<-__________--- Read 3

Read 1 and 2 both have a mate in either exon, so that's good evidence that they come from a single transcript. You can also have reads like Read 3, which spans two exon-exon junctions, again good evidence. There are some problems, however, if things get more complicated: if the cumulative length of your exons is greater than the fragment length of your short reads, or you have many exons in a gene, you will never get a single read pair that confirms the existence of some transcripts. The same goes for the opposite. If you find a read pair on a single exon, that is not evidence that this single exon is a transcript in and of itself.
Here's where short read assemblers take different approaches. Some make educated guesses about the existence of some transcripts, e.g. by utilizing expression counts in an EM or dynamic programming approach, some just enumerate possibilities.
And that's where IsoSeq comes in. A single, long read from a single transcript helps a lot to add evidence and make sure your transcripts actually represent the transcripts that exist in your organism. So, often, your IsoSeq should be fewer transcripts than your short read assembly.
In essence, what you should do, is take the IsoSeq and curate those as a set of high-confidence transcripts, i.e. those that you have observed from a single read source. That's not to say that some transcripts that exist only in the short read assembly aren't real, just that you have less definitive evidence for them.
If you are mostly interested in the gene level, that doesn't make a huge difference, but once you need to look into transcript level usage/expression etc it helps a lot to be able to eliminate low confidence transcripts.
Use tools better suited for transcriptome assembly evaluation:
Detonate, TransRate and (Busco) and compare your assemblies using those.


Find similar proteins with similar tertiary structure?

As a general question on the board, more detail is often good, because of the answers can be inaccurate due to misunderstanding of the biological context.

There are several ways to do this and at least 3 that I'm aware of.

Blast - sequence similarity searching the PDB database
Contact map identity, I use this approach, but its complicated - and its precision or utility is context specific, i.e. it could be loads of work and miss the vital detail you need.
Direct mapping one protein structure against another, this in particular is true for specific domains, or sites within the protein and the result is angstroms. Its not difficult and I assume this can be scaled to the whole protein.

Point 1 "Blackboxes"

By far the easiest is option one and I suspect I-Tasser (here) use this principle before instigating "threading". This package is considered "good" (others appear to consider it the best of its type), I believe it has done well in CASP tests (Critical Assessment of Techniques for Protein Structure Prediction). So the simple answer is dump your sequence in there an watch the output. Pyre is something similar. These are interesting tools but certainly not definitive, really there a starting point.

Point 1 "Formal Blasting PDB"

This is the easiest and most mature approach for your question, without knowing the details. I've run a test to ask the following question.

Which protein structure available on PBD is most similar to Japanese encephalitis virus (JEV) based on protein homology?

Step A
I've no specific interest in the virus BTW. The most robust Blast engine IMO is NCBI here . I dump a JEV sequence into the box, and click three important options (otherwise it doesn't work).

Database = protein databank (PDB), i.e. structural data
Organism "Japanese encephalitis virus (taxid:11072)": the taxid will come up automatically (if not just proceed as normal)
Next to the Organism box I click "Exclude", otherwise this whole approach doesn't work
Obviously I select "blastp", otherwise it definitely will not work

The output is/was here (but is now inactive). Summarising what it presented, E is so high its "0" (>2e-139) protein ID "77.59%" and just states "envelope protein" Chain A. I hit the link and the protein that comes up has a name which ends in "6A0P_A". This is the PDB reference and Percentage ID and PDB reference are the only two bits of info you need. The PDB reference needs the "_A" removing, thus its 6A0P. The _A just refers to chain A of a dimer.

Step B
You then go back to the Web site you mentioned and input the PDB ID you stated into the search bar. This results in this quite interesting page here. The structure is Usutu virus, which surprised me.

Point 2 . This is much harder and not only involves specific Python code, but a post unsupervised learning analysis and results in an average similarity of structure over the whole protein. There are a number of published algorithms for generating a contact map, or its easy to write your own. Its more advanced, its complicated and unless you've seen loads of them they not immediately understandable. Contact maps angle at a very specialist approach of predicting and understanding protein structure. It is not for the faint hearted and some biological details will be omitted, which may or may not be important for your purpose.

Point 3 . I don't have research grade expertise here, so I forget the exact button to click, but its quite easy via PyMol and other graphical viewers to perform spatial similarity for specific epitopes and this is something you need to learn within your sphere. Basically you will load two proteins into the image viewer and then they will be a "similarity aligning" function and you zoom in on your output to generate an answer in angstroms. You will need to learn this approach for the later parts of your project.

Summary, The orthodox way would be to perform point 1, select several targets and then examine them in detail via point 3. I definitely recommend direct Blast if you go down this route, because it is the most tractable and avoids potential unknown subsampling.  I personally would not recommend Point 2 in a general context because you have to consider an unsupervised learning approach towards processing a heap of contact maps. Ultimately, which ever approach you take you can never avoid Point 3 (which someone else does, in my case).

Note This may or may not be obvious to your but just because two sequence have the closest homology doesn't mean they will have the closest structural similarity.


Finding the location and unit length of repetitive sequences within a long read

It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?


Bioinformatics is a big and very heterogeneous field with a lot of variation, so it's hard to recommend something that covers everything. That said, I think the best you'll do is to build familiarity with some of the common tools with some canned analyses.

One possible option would be some of the scRNA-seq workflows from the Pachter group, they seem to be fairly conscientious about documentation. Here is one such canned analysis, here is the paper, here is a blog post about it that suggests next steps and discusses at a high level.

scRNA-seq is just what's been trendy for a few years, next year it will be something different. But it at least gets you familiar with some of the file formats and concepts such as alignment, etc.

Hope that helps.


In short the answer is no, however the 'car-jumping method' described here it is possible, but it is about the enzymic MMLV reverse transcription in response to the capping structure. I assume that an alternative reverse transcriptase doesn't do 'car jumping' therefore any sequencing method will pick this up, but the difference between MMLV verse the heat-tolerant style reverse transcriptases. I would assume the secondary structure associated with viral capping prevents a low temperature reverse transcriptase, i.e. MMLV from correctly reproducing the vRNA into cDNA.
It is important to note that a cap, isn't just a virus thing, its mRNA eukaryotes per se. The cap is a 7-methylguanosine (m7G) linked to the 5′ end of the transcript by a 5′-5′ triphosphate bridge [m7G(5′)ppp(5′)N- or m7G(5′)ppp(5′)Nm in which Nm is a 2′-O methylated nucleoside].


What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?

There are also third party free and open source basecallers that haven't been developed by Oxford Nanopore.

Of particular note is Chiron, which gave the best uncorrected assembly identity among the base callers that Ryan Wick tested. It's slower than Albacore, but appears to be more customisable, and could theoretically be used to model any DNA feature that is detectable in the electrical signal, rather than just the A/C/G/T base sequence.

There's a paper about Chiron on BiorXiv.


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

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

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


Download protein structure from gene name

There may also be a misunderstanding, from your previous post. To reiterate, I don't use Ensembl, which is the source of your gene names.

On Genbank there is an entire "Structural" database for every structure that is on PDB. When Genbank hold the protein sequence that was used to generate the structure Genbank only use PDB codes. For example "5LBV_A" is both a PDB code and a Genbank code. So this means PDB entry 5LBV chain A. Thus its easy to download the crystal because you simply ask PDB for 5LBV. To prove it,

efetch -db protein -format gpc -id "5LBV_A" | xtract -pattern INSDSeq -group INSDSeq -element INSDSeq_sequence

OUT
    ircigvsnrdfvegmsggtwvdvvlehggcvtvmaqdkptvdielvtttvsnmaevrsycyeasisdmasdsrcptqgeayldkqsdtqyvckrtlvdrgwgngcglfgkgslvtcakfacskkmtgksiqpenleyrimlsvhgsqhsgmivndtghetdenrakveitpnspraeatlggfgslgldceprtgldfsdlyyltmnnkhwlvhkewfhdiplpwhagadtgtphwnnkealvefkdahakrqtvvvlgsqegavhtalagaleaemdgakgrlssghlkcrlkmdklrlkgvsyslctaaftftkipaetlhgtvtvevqyagtdgpckvpaqmavdmqtltpvgrlitanpvitestenskmmleldppfgdsyivigvgekkithhwhrsgstiggpfeddddkagwshpqfekgggsgggsgggswshpqfek

What I cannot easily do, is go from a generic NCBI Genbank protein code to PDB code.

However, PDB uses Uniprot codes for proteins, so you may be able to use a Uniprot code to get hold of the PDB code, or Ensembl codes may, like Genbank, use PDB accessions.

So in summary perhaps you need to give more detail?


Test to determine if two genes/exons share the same evolutionary histories?

Since there hasn't been any other answers, I'd thought to share my approach in case someone is having the same issue.

Gauging from the comments of tendon and NatWH, it seems using the SH-test valid stategy to tackle this problem.

Having done some research, I could only find one functional software package, the R-package phangorn which implements both, the SH- as well as SOWH-test and thought to share my R-workflow:

Reading in multiple-alignment (MLN). Here: amino-acids (AA)

 library(phangorn)
 mln <- read.phyDat("my_mln.fasta", format = "fasta", type = "AA")

Partitioning MLN into two parts (part A and B) given a predefined breakpoint (here 100)

breakpoint <- 100
 # partA: 1:breakpoing
 ( mln_partA <- subset(mln, select = 1:breakpoint, site.pattern = FALSE) )
 # 58 sequences with 100 character and 45 different site patterns.
 # The states are a r n d c q e g h i l k m f p s t w y v

 ## partB: breakpoint:end
 ( mln_partB <- subset(mln, select = (breakpoint + 1):length(attr(mln, "index")), site.pattern = FALSE) )
 # 58 sequences with 1702 character and 720 different site patterns.

Constructing tree(s) separately from part A and B using Neighbourhood Joining (NJ) method

 ## Note, other tree building algorithms might be used
 ( treeNJ_partA <- NJ(dist.ml(mln_partA)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.
 ( treeNJ_partB <- NJ(dist.ml(mln_partB)) )
 # Phylogenetic tree with 58 tips and 56 internal nodes.
 # Unrooted; includes branch lengths.

Using topology of previously constructed trees and compute optimal likelihood based on partA data

 ( treeNJ_partA_fit1 <- optim.pml(pml(treeNJ_partA, mln_partA)) )
 # optimize edge weights:  -1402.126 --> -1379.043
 #
 #  loglikelihood: -1379.041
 #
 # unconstrained loglikelihood: -302.6292
 ( treeNJ_partB_fit1 <- optim.pml(pml(treeNJ_partB, mln_partA)) )
 # optimize edge weights:  -1633.646 --> -1530.745
 #
 #  loglikelihood: -1530.682
 #
 # unconstrained loglikelihood: -302.6292

Performing SH-test for for both trees

 ( shPart1_t12 <- SH.test(treeNJ_partA_fit1, treeNJ_partB_fit1, data = mln_partA, B = 1e5) )
 #      Trees      ln L Diff ln L p-value
 # [1,]     1 -1379.041     0.000 0.48461
 # [2,]     2 -1530.682   151.641 0.00000

As can be seen in this example, the tree constructed for partB doesn't fit partA data as much as the tree constructed from partA.


Can FASTA files have nucleotide and protein sequences within them; or must they only have 1 type?

While there's nothing stopping anyone from doing that with the FASTA format (after all, it's just a text file with '>' defining header lines), I don't know of any software that would support such a file structure. At best, it would interpret the nucleotide sequences as protein sequences (A/C/G/T are all valid 1-letter protein codes).

A better question to ask would be "Does ultra-cool-bioinformatics-tool X support combined nucleotide and protein sequences in the same FASTA file?" In which case the answer would most likely be, "No."


How to correctly call a VCF file using damaged DNA? (IonTorrent & FFPE)

Sites that were deaminated are often called abasic sites.

Q1: How to obtain DNA damage metrics?

You might want to check out mapDamage - this software is used to perform these sorts of analyses and produce plots of types of DNA damage along sequencing reads.

Q2: How to correct for abasic sites and other DNA damage errors?

You can use your choice of kmer-based or alignment based error correction software. The errors will be random and if the sequencing depth is sufficient the erroneous kmers at the ends of reads will be corrected. However, depending on what you are trying to do downstream with your reads you may be better served by not correcting the reads, then dealing with transitions/damage later.


Toptable error, wont recognize condition

The error is telling you that there isn't any Infected column on the model.matrix.
Check the column names of mm. Also not sure that having two ~ will work well.


How to calculate lipophilicity of proteins?

Oligopeptide

The lipophilicity of an oligopeptide or a small molecule can be calculated with RDKit, the classic compChem python toolkit. The measure is logP or partition coefficient, namely the magnitude of the ratio of the concentration in octanol over than in water, so the higher it is, the more it partitions in the octanol phase. However, with a compound of several dalton it may be good to see what the contribution of each atom to this score are. This measure is the Crippen contribution to logP, which can be handily plotted on the molecule. Note that these are may be at first glance look similar to Gastaiger partial changes, but are not.

Okay, first the imports

from rdkit import Chem
from rdkit.Chem import AllChem, Descriptors
from rdkit.Chem.Draw import SimilarityMaps

Now let's make the molecules. Let's start with something we can simply use a SMILES string of the web: Bortezomib.

name = 'Bortezomib'
mol = Chem.MolFromSmiles('O=C(N[C@H](C(=O)N[C@H](B(O)O)CC(C)C)Cc1ccccc1)c2nccnc2')

Now we can do all the relevant calculations on the mol object:

print('logP', Descriptors.MolLogP(mol))
AllChem.Compute2DCoords(mol)
contribs = Chem.rdMolDescriptors._CalcCrippenContribs(mol)
fig = SimilarityMaps.GetSimilarityMapFromWeights(mol, [x for x,y in contribs], colorMap='BuPu', contourLines=10)
fig.savefig(name+'_crippen.png', bbox_inches='tight')

LogP=0.3

colorMap='BuPu' is the matplotlib colorscheme. I have no idea which work best, I have spent way too much time trying and I do not know. But with BuPu, white is hydrophilic, purple is hydrophobic.

Say the oligopeptide does not have a SMILES or we don't want to do it that way. Let's take β-Casomorphins 1–3:

name = 'β-Casomorphins 1–3'
mol = Chem.MolFromFASTA('YPF')
mol = AllChem.ReplaceSubstructs(mol,
                                Chem.MolFromSmiles('N'),
                                Chem.MolFromSmiles('[NH][OH]'),
                                replacementConnectionPoint=0)[0]

The result with the last block from before gives:

logP=1

Proteins

With protein, it is a lot harder. Membrane affinity of protein is dependent on their structure. In PyMOL you can use the Adaptive Poisson-Boltzmann Solver (APBS) to show the surface by charge. For transmembrane protein there is clear ring of hydrophobicity —CHARMM-GUI membrane builder and a variety of tools for example can guess where it is. For a membrane associated protein only one side will be hydrophobic with a tell-tell charged ring around it (e.g. DHR1 domain) and there are no hard and tested methods for that.


Running snakemake with multiple conda environments without rebuilding environments every time

Each YAML file will always create its own environment. The packages in each aren't downloaded and installed each time, if they're already installed elsewhere they're just hard-linked in. You'll want to use both --use-conda and --conda-prefix. The latter will point to the base directory in which all conda environments will be stored. This has the benefit that running the same workflow in different places won't result in a new env being created (by default snakemake will do so inside a .snakemake directory otherwise).


Any user friendly way to find rare mutations in whole genome raw?

You need to run a SNP finding pipeline, filter for SNPs in those genes and then run a SNP annotation pipeline. (You could annotate before you filter but it will take longer).

SNP finding usually has three steps, 1) read alignment to reference, 2) SNP calling and 3) SNP filtering (by quality).

You then bring in gene annotation data to filter SNPs to just those genes.

You then run SNP annotation using tools like VEP, SIFT, PolyPhen, etc.

People use tools like GATK to do this.

Hopefully that's enough information to get you started.


Adjusting phenotypes by regressing out covariates

Instead of the "observed phenotype", it's asking for a "corrected phenotype", which is the residuals from pheno ~ covariates. In R, one would get that as follows:

fit = lm(phenotype ~ covariates, data)
corrected_phenotype = resid(fit)

Of course covariates is just a place hold for something like age + gender or something along those lines.


Looking for a taxonomic dataset that links tree genera (or species) to families

The NCBI taxonomy is a quite extensive resource: https://www.ncbi.nlm.nih.gov/taxonomy

I would thus expect there is a good chance that it covers all of the taxa that you are interested in.  It officially claims that it is not an authoritative source, but it's still pretty darned good.


How to read the "SEQRES" section from a PDB file, using R

I have not used Rpdb but I have used, and would recommend, bio3d.

To read the PDB file 1A4G:

library(bio3d)
pdb1a4g <- read.pdb("1a4g")

To get the SEQRES information, there are two options. Get it from the pdb object directly:

seq1a4g <- pdb1a4g$seqres

Or call pdbseq:

seq1a4g <- pdbseq(pdb1a4g)

Using the second method, the 17th position (note this is not the same as the 17th position of the original protein sequence):

seq1a4g[17]
 92
"Q"


How can I subset a reference based on only the first chromosome?

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


Where can I visualize Ramachandran plot?

Here is a page with several suggested methods. I agree that Rampage looks the best of the options, but it seems that it provides some suggested workflows for drawing it in R so that it is very customizable and you can probably get it to show any features you want if you're willing to customize a little.

Here is a Python tool for drawing Ramachandran plots.

There are additionally several other servers/tools.


Display a score or graph to visualize the degree of conservation of each residue from alignment data

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

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


How to quantile normalization on RNA seq counts


gene-gene correlation from two different Tissues


Quantification of EST/cDNA based retained intron transcripts

You can do full-length whole-transcriptome sequencing using either the direct RNA kit or a cDNA kit (with or without PCR) on a MinION. This will allow you to quantify isoforms directly, rather than estimating by looking at read pairs.

However, it'd be a bit expensive to do it that way if you're looking at a single transcript (and aren't barcoding to multiplex samples on the same run).

If you've already done Illumina whole-transcriptome sequencing (e.g. strand-specific paired-end sequencing), tools like Kallisto, RSEM, or Salmon can be used to estimate isoform expression based on the mapped reads. These tools use the information about transcript placement and/or similarity to alter the observed read counts in an attempt to better represent reality.

Unfortunately, if you want these tools to quantify intronic sequence (i.e. the $2^4$ possibilities you have suggested), then those sequences need to be included in the sequences provided. They won't estimate frequencies for non-existent sequence.

As an alternative to including every possible gene model, something like a SuperTranscript model might work. I'm not familiar with Alicia Oshlack's method, but it at least gives the appearance that it might reduce the complexity of searching for multiple unannotated isoforms within a gene.


samtools mpileup empty when filtering out flags

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:

...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):

...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l


Reference genome for allele specific expression

As you said, this is a well-known problem for ASE. In Degner et al. 2009, they state that, "Perhaps surprisingly, we found that masking known SNPs does little to eliminate inherent biases in read-mapping."

As such, substituting N bases is not likely to completely solve your issue. Depending on the aligner, these N bases may also negatively affect your ability to map reads.

Of your two suggestions, I believe HISAT2 will be the better approach. There are also more intensive methods you can use such as empirical Bayesian methods or more traditional variant phasing (ex. AlleleSeq).


Mix globbing and wildcards when specifying rule input

A nice solution comes from a colleague of mine, as seen in this Snakemake workflow. The trick is to access the wildcards programmatically using an anonymous (lambda) function. In my example above, it would be implemented as follows.

rule preprocess:
    input: lambda wildcards: glob('sample-{samp}/*.fastq'.format(samp=wildcards.samp))
    output: '{samp}-clean.fastq'
    shell: 'mypreproccmd {input} > {output}'


Estimate the length of poly-A tails from randomly-primed RNAseq data

Can’t be done. If you already sequenced then I’m afraid the money is wasted (unless of course the data is good for something else).

The standard Illumina basecaller doesn’t deal well with homopolymers and will consequently call incorrect poly(A) lengths.

Amongst others, Narry Kim’s lab developed a dedicated method — TAIL-seq — to solve this problem. It’s a protocol for library prep using pulldown of biotin tagged tails (rather than random priming), plus dedicated basecalling software that uses spike-ins for calibration.

Be warned though: the protocol is reportedly technically challenging, and in order to use it you need an Illumina machine that gives you access to the raw image files before basecalling. Up-to-date Illumina machines won’t give you these files, full stop. So you’ll need to find a machine that runs an outdated firmware.


Download complete Phylogenetic Tree as JSON

NCBI maintains an SQL dump of their taxonomic database
You can see what that contains here
For smaller requests (depending on what you need), you might also consider the ENSEMBL REST API - although I don't know of a way to retrieve the complete tree (shoot them an email): https://rest.ensembl.org/taxonomy/classification/53399?content-type=application/json


NCBI xml-format Printing nodeValues from tags with the same names

use a XSLT processor (http://php.net/manual/en/class.xsltprocessor.php) with the following stylesheet.

<?xml version='1.0' encoding="UTF-8"?>
<xsl:stylesheet xmlns:xsl='http://www.w3.org/1999/XSL/Transform' version='1.0'>
    <xsl:output method="html"/>

    <xsl:template match="/">

            <xsl:apply-templates select="//Gene-commentary_heading[ text() = 'GeneOntology']"/>

    </xsl:template>

    <xsl:template match="Gene-commentary_heading">

<xsl:value-of select="name(.)"/>

<xsl:value-of select="."/>

        <xsl:apply-templates select="..//Other-source_anchor"/>
    </xsl:template>

    <xsl:template match="Other-source_anchor">

<xsl:value-of select="name(.)"/>

<xsl:value-of select="."/>

    </xsl:template>

</xsl:stylesheet>


Is there a way to import tables in PubMed Central to MySQL or other Relational databases?

If you can use R, consider using web scraping.

See here : http://bradleyboehmke.github.io/2015/12/scraping-html-tables.html


Metapath2vec on Drug-ADR Heterogeneous Graph

Clearly the additional references do put Metapath2vec into good use. Its a many to many ML relationship. Most catagorization/classification particularly in deep learning is many  to one. So thats the good thing about the technique you  are using, you've 25 catagories, which is quite large for machine learning(ML) and large for deep learning. The technique appears to be half way between supervised and unsupervised learning. The random walk establishes a zone (circle) of catagorisation.
NLP Background
The situation where it would work in NLP (natural language processing) is for example a spell checker, where you have one correct word and loads of possible typos, some of which you not observed in training. So what will happen is the random walk will pick up the local vacinity around the 'correct spelling' by the vectorization (converting a word into numbers) of the incorrect spelling will land it close to 'true spelling' in a 2D vectorization plot. The training following the random walk says  'yep you belong with me'.
The papers given talk about being able to identify the context of a word within a sentence which is a very NLP thing, one vectorisation being 'bags of words'. So here the sentence represents one of many categories of information, i.e. each is the central node within the graph, the vectorisation of the sentence places of meaning closer to one 'hub' than another and again the local vacinity is identified by the random walk. In training this 'random walk vacinity' is identified, so it will algorithmically recreate the 'random walk locality' so that sentences that have not been trained can have their meaning identified.
The central theory of deep learning, which not quite true of other ML, is any algorithm can be recreated through training. So like it will sort of recreate a 'random walk zone'
Difference with drug discovery
The papers you have map drugs to genes, but it means you have to have a "true drug to gene" relationship. The problem I see here is what is 'true' for a drug-protein interaction? With a spell-checker there is one 'true spelling', with a sentence there is pretty much one 'type of meaning', like a movie was great, average, terrible based on a ratings score and stuff like gender classification. If its not 'true' the random walk will get it wrong. So its a big assumption. If you have a circle and pick a point on the edge of the circle then try and draw a circle around the original circle, you can't unless the circle is so large it risks including stuff that is not part of the original circle. The problem the papers is how do they know what is the 'centre' of the circle? Its an a priori assumption and you need really good biological rationale.
With your stuff I agree its easier because adverse drug reaction, some drugs could be fatal or toxic (traditional cancer drugs or parasite drugs), so you'd pick the worst case adverse reaction. That then centres the definition around which everything else is relative. If however your 'true catagory' turned out to be an outlier its not going to work, because the random walk only works if its starts bang in the centre of the 'catagory'. Random walks form spheres, in case you are not aware. I can see this working as an advanced 'unsupervised learning' technique.
One very real advantage is if a drug doesn't have an adverse reaction it MIGHT (I'm guessing) not give a result, in which case you are good to go - that a good drug. If that is what happened that would be amazing because its not how ML, or particularly DL works. Traditional ML/DL will force a classification even if it doesn't fit. There are methods to accept some sort of "half-classification" but I don't know how to do it. This is one of the reaons you need lots of things to train.
I am assuming that vectorisation of your drug is fine because there is a huge history of this in drug-discovery going right back to the origins of the field.
Summary
So yeah I get the idea, I see the flaws and the advantages. The final issue is that most ML requires pretty large data sets and DL require huge data sets and when you've got 25 catagories. For DL approaches the data set size might need to be huge unless you can do 'transfer learning' (backing your training onto a pre-trained model). If your method resulted in a zero classification i.e. "doesn't fit anything I've got" - that would be cool and that is the key advantage I can see in the 'random walk' approach which is not available using traditional methods.
As long as you are aware, you do not have to use this algorithm to get a classification. This final point what appears missing in the papers is using a 'traditional' method by comparison. The ML before deep learning, stuff like random forests, routinely takes loads of catagories. If you are using a non-deep learning ML you would traditional use several different approaches to justify your result.

My answer is a bit rambly but I think you'll get the gist. Peronally I would complement with a standard ML, but all methods have their flaws. Deep learning (DL) has a big history in drug discovery (alot of work was done via Tensorflow), but I'm not sure you've data set size for this.


KeyError when getting features from a genbank file with biopython with some accessions but not others

Good question.
feature.qualifiers is a dict, if the dict doesn't have that key it will throw a KeyError. The way your code works it expects feature.qualifiers to have 'gene' as a key every time feature.type == 'CDS'. However, it is possible to have feature.type == 'CDS' without a gene entry and in the case of AP019703.1 this comes up before you reach the gene acrA.
As simeple fix is to simply use get instead:
feature.qualifiers.get("gene",[])

Now returns a empty list if there is not key "gene", which also satisfies the 'name in' part of your if statement.
If you get stuck on stuff like this use more print statements, break up your if statements and look at the whole file (example below) to understand the structure.
gi_genome = 'AP019703.1'
name = 'acrA'
Entrez.email = 'email'
handle = Entrez.efetch(db="nucleotide", id=gi_genome, rettype="gbwithparts", retmode="text")
record = SeqIO.read(handle, "gb")
handle.close()
with open('AP019703.gb' ,'w') as fout:
    SeqIO.write(record, fout, 'genbank')


Is there a safe catch-all adapter sequence for trimming?

You're best off just using fastp or Trim Galore!, both of which will determine the adapter sequence for you. Trim Galore! uses a built-in list of known sequences for this, whereas fastp uses read overlap.


Are there free cloud computing platforms for biology projects?

If you're getting started with a project and need some free computational resources you should look at https://galaxyproject.org/. There are some free servers that you may be able to use. However, Galaxy is designed to have pre-fabricated pipelines and a GUI to avoid having to deal with command-line tools. Therefore any sort of custom analysis will have to go through the galaxy API.

Otherwise, most of the free computational resources that were available several years ago like iPlantCollaborative are no longer available. Running servers costs money, after all, and the funding for such projects often runs out after a few years.

You may consider asking researchers at your local universities to collaborate and provide you with access to their servers.


It is utterly arbitrary and depends on your resources.
We can all agree your cutoff would be a lot lower during the ordering step than at followup in silico analyses, but even then if you were ordering the top N at \$100-250 each, should you spend \$1,000 or \$10,000?
In the in silico step you have real time and CPU-time to factor in making this call even blurrier.
But there are some things to put those scores into context:

The units are kcal/mol, i.e. these are predicted ∆∆G of binding. Here is an explanation of what does a kcal/mol mean for docking — 1 kcal/mol is a hydrogen bond or the average collision energy of water at 37°C.
If you have a crystal structure with a known ligand, re-docking and cross-docking that would give you a benchmark of what you should expect ——keeping my mind that natural ligand are horrid and don't obey Lipinski's rule. The natural ligand will wander and to assess how good the docking went RMSD is used, but for this I mean looking at its GlideScore as if it were any other compound.
Rescore with a different tool —Gold etc.

Ligand efficiency
An important thing to remember is that ∆∆G is size dependent —a huge peptidomimetic will score more than a smaller molecule that hits the hotspots perfectly. Therefore, ligand efficiency (LE), ∆∆G over heavy atoms (reviewed in detail here) is a good metric to prevent huge molecules dominating spuriously. If your library has a variety of different sizes, it is best to sort by LE than ∆∆G alone. In fact, if you get a few small hits you can merge/link them or expand on them, whereas if you get nothing because you gambled on larger hits you have to start from square one. I would argue that this is true for repurposing screens and the low success of repurposing VLS despite the enthusiasm is due to the neglect of this.


What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)?

If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.
What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.
There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.


How can I create my own GO association file (gaf)?

Here's a Perl script that can do this:

#!/usr/bin/env perl
use strict;
use warnings;

## Change this to whatever taxon you are working with
my $taxon = 'taxon:1000';
chomp(my $date = `date +%Y%M%d`);

my (%aspect, %gos);
## Read the GO.terms_and_ids file to get the aspect (sub ontology)
## of each GO term.
open(my $fh, $ARGV[0]) or die "Need a GO.terms_and_ids file as 1st arg: $!\n";
while (<$fh>) {
    next if /^!/;
    chomp;
    my @fields = split(/\t/);
    ## $aspect{GO:0000001} = 'P'
    $aspect{$fields[0]} = $fields[2];
}
close($fh);

## Read the list of gene annotations
open($fh, $ARGV[1]) or die "Need a list of gene annotattions as 2nd arg: $!\n";
while (<$fh>) {
    chomp;
    my ($gene, @terms) = split(/[\s,]+/);
    ## $gos{geneA} = (go1, go2 ... goN)
    $gos{$gene} = [ @terms ];
}
close($fh);

foreach my $gene (keys(%gos)) {
    foreach my $term (@{$gos{$gene}}) {
        ## Warn and skip if there is no aspect for this term
        if (!$aspect{$term}) {
            print STDERR "Unknown GO term ($term) for gene $gene\n";
            next;
        }
        ## Build a pseudo GAF line
        my @out = ('DB', $gene, $gene, ' ', $term, 'PMID:foo', 'TAS', ' ', $aspect{$term},
                             $gene, ' ', 'protein', $taxon, $date, 'DB', ' ', ' ');
        print join("\t", @out). "\n";
    }
}

Make it executable and run it with the GO.terms_and_ids file as the 1st argument and the list of gene annotations as the second. Using the current GO.terms_and_ids and the example annotations in the question, I get:

$ foo.pl GO.terms_and_ids file.gos
DB  geneD   geneD       GO:0005634  PMID:foo    TAS     C   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0003677  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030154  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
Unknown GO term (GO:0006350) for gene geneD
DB  geneD   geneD       GO:0006355  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0007275  PMID:foo    TAS     P   geneD       protein taxon:1000  20170308    DB
DB  geneD   geneD       GO:0030528  PMID:foo    TAS     F   geneD       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016020  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005524  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006468  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0005737  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0004674  PMID:foo    TAS     F   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0006914  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0016021  PMID:foo    TAS     C   geneB       protein taxon:1000  20170308    DB
DB  geneB   geneB       GO:0015031  PMID:foo    TAS     P   geneB       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0006950  PMID:foo    TAS     P   geneA       protein taxon:1000  20170308    DB
DB  geneA   geneA       GO:0005737  PMID:foo    TAS     C   geneA       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003779  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0006941  PMID:foo    TAS     P   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005524  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0003774  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005516  PMID:foo    TAS     F   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005737  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB
DB  geneC   geneC       GO:0005863  PMID:foo    TAS     C   geneC       protein taxon:1000  20170308    DB

Note that this is very much a pseudo-GAF file since most of the fields apart from the gene name, GO term and sub-ontology are fake. It should still work for what you need, however.


How to use C++ htslib to read VCF contig name and size?

#include "htslib/vcf.h"

int main(int argc, char *argv[])
{
    htsFile *fp;
    bcf_hdr_t *hdr;
    bcf_idpair_t *ctg;
    int i;
    if (argc == 1) {
        fprintf(stderr, "Usage: print-ctg <in.vcf>\n");
        return 1;
    }
    fp = vcf_open(argv[1], "r");
    hdr = vcf_hdr_read(fp);
    ctg = hdr->id[BCF_DT_CTG];
    for (i = 0; i < hdr->n[BCF_DT_CTG]; ++i)
        printf("%s\t%d\n", ctg[i].key, ctg[i].val->info[0]);
    bcf_hdr_destroy(hdr);
    vcf_close(fp);
    return 0;
}

On stability: this use has been in htslib forever. In general, functions/structs/variables in the public headers are meant to be stable. However, there is no guarantee that future versions will always keep the same APIs.


How to remove degenerate sequence (N4) at the 5' and 3' of the reads?

From the link you posted it looks as though they used UMI tools as per the section which says
'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'.  The link above has instructions for installation and use.
From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read.
Looking up the protocol for the nextflex kit used (See Link) confirms this:

  Data Analysis
  The 3' and 5' adapters included in this kit both contain 4 random bases that will appear immediately 5' and 3' to the insert in sequencing data. The presence of these random bases should be considered when choosing an alignment strategy. When using “end-to-end” alignment, we
  recommend processing data in the following manner:
  1. Clip the 3' adapter sequence (TGGAATTCTCGGGTGCCAAGG).
  2. Trim the first and last 4 bases from the adapter-clipped reads.
  3. Perform alignments as normal.

It seems that N4 is a 4 base random barcode at both the 5'& 3' end of the read - these are cut from the read and this metadata is added to the read name so that every read can be uniquely identified.

Edit:

You can install UMI tools via:

conda install -c bioconda -c conda-forge umi_tools

or

pip install umi_tools

Usage instructions can be found at:

https://umi-tools.readthedocs.io/en/latest/reference/extract.html

  Usage:
  For single ended reads, the following reads from stdin and outputs to stdout:

  umi_tools extract --extract-method=string --bc-pattern=[PATTERN] -L extract.log [OPTIONS]

  For paired end reads, the following reads end one from stdin and end two from FASTQIN and outputs end one to stdin and end two to FASTQOUT:

  umi_tools extract --extract-method=string --bc-pattern=[PATTERN] --bc-pattern2=[PATTERN] --read2-in=[FASTQIN] --read2-out=[FASTQOUT] -L extract.log [OPTIONS]


Omics data: How to interpret heatmap and dendrogram output?

The dendrogram summarize the information of a group of values and sort them according to the similarity they have. It can be applied to both, samples and features.

The dendrogram allows to visualize features that are more similar together, usually revealing patterns that wouldn't have been seen otherwise.

In an article usually it is used something like: "we can see that X and Y are clustered together revealing that ... (they have something in common), while Z and A are clustered far"


Can we use GTEx data as control data for TCGA data?

Definately not directly. Probably not even if used indirectly. Different projects, done at different times, with different protocols, different machines, different library preps etc.

We have tested combining GTEx and TCGA data for several different purposes where one might expect the analyses' internal controls would have allowed compatibilty, but even then we have always found a persistent batch effect which we were unable to remove.


Determining the quality of a publicly available genome assembly

It is not so easy to determine the quality of an assembly - but sequence length, number of contigs, N50 and L50 can give you some guidance.

How much does the sequence length deviate from the expected length (or how much variation is among assemblies of the same species). Same for number of contigs/scaffolds. The sequence length, and number of contigs together with N50 and L50 give you an idea about the fragmentation. Very small N50 compared to expected genome size and chromosomes is a likely indication of high fragmentation, same for high L50 compared to expected number of chromosomes.

However, these number can be misleading. High N50 and low L50 can occur witho some very long contigs, so you may want to compare N75/L75 or better the Quast statistics


Heteroplasmy at any site (as long as it's a reliable variant) suggests that there is mitochondrial heteroplasmy, it doesn't matter where it is. In that case, the "mitochondrial heteroplasmy" statistic would simply be the maximum value at any [reliable] location.

Due to the low amount of variation in the mitochondrial genome, and because it's a haploid genome, I don't think that a calculation like "average heteroplasmy" makes much sense for mitochondria.


Inferring a phylogenetic tree from BLASTn

If you look into the get_seqrecs() function, you will see that this function is not extracting hsps to use in the following steps but the hsps are used in order to filter out low score BLAST alignments. The function accepts a blast alignment result and a threshold, and if the hsp score is above the supplied threshold the whole sequence (not just the hsp) would be the output. Your intuition is correct, the whole sequence should go in the MSA step.


Generating random BED intervals given constraints

You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6).

To help with this, you could use fetchChromSizes from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest.

Given a chromosome and a maximumEnd that falls within its bounds, uniformly sample from [0, maximumEnd-minimumStart-maxLength), adding minimumStart to determine start position and maxLength as length to get the stop position.

To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into bedmap --echo-map-size operations on repeatmasked regions (4) and bedops -n 1 on "regions-to-avoid" (5) operations to do rejection sampling.

Perhaps something like:

$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1<threshold)' \
    | cut -f2- \
    | bedops -n 1 - regionsToAvoid.bed \
    | head -${SAMPLES} \
    > qualifyingCandidates.bed

This will be relatively efficient, because once the head process at the end writes the number of SAMPLES to qualifyingCandidates.bed, it will send a SIGINT signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.

Hope this helps!


Inflated p-values in quantitative trait analysis

Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.

If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.

Some random thoughts:

What does the distribution of your methylation data look like after normalization?
Your observed distribution is basically never on the expected line, which argues strongly for confounding.
How are you adjusting for non-independence of methylation types? Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model, or at least fit 10-20 principal components to account for population structure.
How are you representing your methylation status? If it's binary you might look at this paper.


plot correlation between several bam files

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

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


How many reads do I need for hybrid assembly

I have never used dbg2olc. However other OLC assemblers managed to do a solid assembly from 20x of corrected long reads. But obviously, the more the better.

In you case of allotetraploid, you need to figure out if the heterozygosity is sufficient to separate haplotypes of not and adjust parameters of the assembler. For instance in bdelloid rotifers the tetraploid genome has an AA'BB' structure, all the assemblies were able to separate A and B, but not A and A' and B and B'. I guess an allotetraploid could have a similar structure.

It would be good to know what you are dealing with and what is your goal in the assembly. Do you want a single haplotype?

-- edit --

You can check the genome structure using short reads by smudgeplot, a tool I made (yes, it's a shameless plug).


Error rate setting in Canu error correction

It's generally a good idea to trust the "official" suggestions. You can also adjust the error rate based on coverage according to parameter reference:

  For low coverage datasets (less than 30X), we recommend increasing
  correctedErrorRate slightly, by 1% or so.

  For high-coverage datasets (more than 60X), we recommend decreasing
  correctedErrorRate slighly, by 1% or so.

  Raising the correctedErrorRate will increase run time. Likewise,
  decreasing correctedErrorRate will decrease run time, at the risk of
  missing overlaps and fracturing the assembly.

That being said, the sequencing quality can vary from between libraries or flow cells. The error rate is not a constant.

According to the developer (full thread):

  if you want to assemble 1D data or have a bad run, it is possible you
  would need to increase the error rate. You could run with a high error
  rate (0.1) and look at the distribution of overlap error rates in the
  unitig step to look for a peak in the distribution. If you have a near
  neighbor you could also map the corrected reads to it and estimate the
  residual error that way.


Weighted Network Analysis Book by Steve Horvath, taking the power of a matrix

You CAN calculate M^b for b being float number, complex number or even another matrix.
As well as you can calculate exp(M), log(M), sin(M) whatever... That is very standard for math people.

Let me give two ways how you can think of it.

1) M^b = exp(b log(M)), so you need to define what is exp and log for matrices.
That can be done using power series: exp(M) = 1 + M + M^2/2! + M^3/3! + ...
log(1+M) = M - M^2/2 + M^3/3 - ...

2) More general way to calculate any function of matrix f(M) , looks like that:
if your matrix is diagonal, then it is true that: f(diag(d_i)) = diag( f(d_i)  ) .
So the question how to reduce any matrix to diagonal one.
That can be done for almost all matrices by conjugation M = C D C^{-1}.
Then f(M) = C f(D) C^{-1}.


What do the first few lines of your lmFit input variables look like (i.e. dat and pheno)? Are the rows numbered, or do they use the actual symbols?

By default, I don't think read.csv expects a file to have row numbers (even though write.csv puts them in by default). Try adding a row.names parameter to the expression matrix read call:

dat <- read.csv("expression.csv", row.names=1)


Here are a couple books I'd recommend:

Dan Gusfield's Algorithms on Strings, Trees, and Sequences is a deep and wide treatment of aligning, searching, and processing strings, trees, and sequences.

Warren Ewens' Statistical Methods in Bioinformatics devotes a chapter to BLAST and the math underneath it.

Edit - Another book that may be useful for some mathematical underpinnings would be Michael S. Waterman's Introduction to Computational Biology.


RapMap: reference transcriptome for simulated reads

Thanks for your interest in RapMap.

At that time we were using flux simulator for simulating read sequence data. We used the genome and gtf file together as an input to flux.

I dug into the scripts and got hold of the gtf link ftp://ftp.ensembl.org/pub/release-80/gtf/homo_sapiens/Homo_sapiens.GRCh38.80.gtf.gz,
and the genome file link ftp://ftp.ensembl.org/pub/release-80/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa.

We used rsem-prepare-reference to create the reference transcriptome file later in order to run RapMap and other transcriptome based quantification tools. I hope this helps.

PS: I am a coauthor of the RapMap paper.


Where can I find the genotype frequencies for the SNP's that are tested at 23andMe, Ancestry and FTDNA/MyHeritage?

gnomad is probably the most comprehensive allele frequency database that exists at the moment.
Depending on how many variants you need to query, you can either manually input them into the search bar, or if you have a lot more, then you can also directly download the vcfs under the downloads tab.
The default for V2 gnomad set is on the gr37 build.


BAM to BigWig without intermediary BedGraph

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

Given single-end file for sample1.

library(GenomicAlignments)
library(rtracklayer)

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

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

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

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


How are Principal Component analyses and Admixture analyses from a genetic alignment different?

Just to add to the great summary by Devon Ryan: admixture analysis tools are much more flexible than PCA (which is just a fixed mathematical operation), so they can be designed to incorporate LD patterns, phase information, different population evolution models etc.

Here's a paper from the creators of fineSTRUCTURE and Chromopainter admixture tools detailing the action of both - it seems they provide some comparisons with PCA as well:
http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002453


pvalue calculation of protein-protein network with permutation test

It is probably necessary to correct for the degree distribution of the network, I am not sure that a tool such as the Mantel test will quite do this (as suggested by user Michael G.). In naturally occurring networks degree distribution is a strong confounder.

Then again I am not totally sure what shape your data is in or what the question is. You mention a background network, but I am not sure what this is or how you want to use it. Possibly your network reconstruction is a subset of your background network? In which case you just define the rule for how your subsample is drawn and then you compare it to permuted networks.

However, here is some fairly crude but serviceable R code that implements a general purpose network rewiring test, using igraph. I am not using a background network as again I'm not sure what it is.

(Python also has igraph so should be similar rewire() function)

compute_measure = function(network) {
  # i don't know what you are measuring, so you have to write this function.
  # this function takes a graph object as input
  # this function returns a single numeric value (whatever you're measuring) as output
}

# "net" is a network of the graph type from igraph

measure = compute_measure(net)

num_perms = 1000

perm_measures = c()

for (i in 1:num_perms) {
  # permute the network by rewiring, rewire 10x as many times as edges
  # to ensure saturation
  perm_net = rewire(net, with = keeping_degseq(niter = gsize(net) * 10))
  perm_measures[i] = compute_measure(perm_net)
}

# one-sided monte carlo p value is proportion of rewirings (plus one) greater than your measure
p_value = (1 + length(which(perm_measures >= measure))) / num_perms

If on the other hand you are just trying to come up with a single permutation of the edges of your network to compare against, you can run rewire() only once to get that permuted network. It depends on what question you are trying to answer, and what null hypothesis you're trying to get a p value for.


How do researchers account for heterozygosity in genome assembly?

Generally a single allele is used. This may be arbitrary or match the major allele frequency of the species. For sex chromosomes, usually both/all are included, so for mammals both X and Y are part of the reference genome.


If you're building a tree of different species, it's probably better to use protein sequences. A common resource is Uniprot.
If you still want DNA sequences you can turn to NCBI Nucleotide.
Some examples that give decent trees: p53, TOP1, HBA.
A common way to construct a phylogenetic tree is to start with a multiple sequence alignment (MSA).


Parse out exon coordinates from bed file for each gene

There are multiple ways to go about it. On the command line you can make a 1 line BED file:

chr1   11868   12227

And then bedtools intersect with it.

In R, you could load your original BED file and use GenomicRanges:

library("GenomicRanges")
bed = read.delim("foo.bed", header=F)  # Rename this
# N.B., BED files use 0-based coordinates, I've switched to 1-based here
exons = GRanges(seqnames=levels(bed$V1),
                ranges=IRanges(start=bed$V2+1, end=bed$V3),
                strand=bed$V4)
# Again, I switched to 1-based coordinates...you don't have to do that
region = GRanges(seqnames="chr1",
                 ranges=IRanges(start=11869, end=12227),
                 strand="+")
theExonsYouWant = subsetByOverlaps(exons, region, type="within")

You'll need to either specify a strand in your region or use ignore.strand=T in subsetByOverlaps(). I hope you don't have any genes that overlap on the same strand, since they'll prove problematic.


scRNASeq expression matrix with decimal values

It looks like you have normalized counts. My guess is that they are log-normalized with some scale factor, probably 10000.

With the whole matrix you should be able to work out the scale factor for each cell. See https://www.nxn.se/valent/2018/10/25/unscaling-scaled-counts-in-scrna-seq-data

They typical format for scRNA-seq data is to have cells in the columns and genes in the rows. So here your rownames are the gene names and colnames are the cell names.


How to share a collection of sbol designs in SynBioHub

There is a balancing act between allowing all edits, which defeats the purpose of having URIs point to static objects, and allowing no edits, which defeats the purpose of providing some updated information about an object.  SynBioHub resolves this issue by allowing a limited amount of edits after publication.  To answer your questions specifically:

Yes, you can share objects by generating a share link OR adding another owner to an object.  Sharing can be done by the original owner of the object.
There are some mutable fields that any owner of the object can change, these include the description, notes, source, and references in the Details section.  You can also add attachments.
Yes, when you publish (Make Public) a new object or collection, you have the option of making it public into an existing collection OR a new collection.  A person with curator privileges on the SynBioHub instance you are working on can perform this operation.
Answered above.


How would we test for significant changes in expression when many of the normalised counts are clustered together?

TLDR- just do the analysis in DESeq2 (it can handle multiple covariates)
If I've interpreted your question correctly there seem to be a few confusions here.
First off, for a multiple linear regression aka ANCOVA, one of the assumptions are the residuals are normally distributed - however this is almost never an issue in practice. The distribution of the data is not so important, you make the assumption that the data has come from a population which is continuous and normally distributed that's why you've chosen to do an ANCOVA. But, RNAseq produces count data - not continuous and bounded at zero. You need to be looking at using a Poisson GLM, or negative binomial GLM as done by DESeq2 (lookup overdispersion)
The normalisation done by DESeq is to account for different library sizes,it doesn't transform the data to make it conform to a normal distribution. You'd need to do some form of data transformation to do that, but don't it's unnecessary.
I'd suggest reading the DESeq vignettes and following their standard pipeline.
Hope that is of some help.


calculateQCMetrics defunct, how to calculate the quality metrics by perCellQCMetrics

Ok I have figured it out and I am sharing this in case anyone can get help. I used
sce <-addPerCellQC(sce)

which calculated the QC for all the altexps also and adds them to the colData.


nextflow `fromPath` not passing all files with *

got it!

I should run:

nextflow run main.nf --input "*bam"

instead of

nextflow run main.nf --input *bam


extract chimeric and multimap reads from bam file

This is how I do it.
samtools view -@5 -f 0x800 -hb /path/sample.bam > sup.bam

Then if you want it as a fasta
samtools view -F 0x1 -hb sup.bam | samtools fasta -F 0x1 - > sup.fa


How to read gene regulatory network edge list files?

You can read about it on the website, the last column is the Log likelihood scores (LLS) between the genes in the first two columns. An LLS of 0 means not better than random link.

On the website they show papers with fuller details about the statistics.


What is a proper way for random subsampling of metagenomic data?

I am not NGS expert, I known ML. Essentially you are performing data augmentation which is essential for ML. What I question is whether a replicate of e6 from a population of e7 would really be sufficient for ML because you need sample sizes (replicates in this case) of >100000.

The answer is really simple you bootstrap the data (sample with replacement), to obtain your replicates. Numpy has a random module that will do this,

p.random.seed(345)
pop = np.random.randint(0,10000000 , size=1000000)

Use pop to select from fastaq. The other way in is to define each fastq sequence within pop then simply,

sample = np.random.choice(pop, size=1000000)

Or a combination of the two. If you can do ML you can do Python coding.


Get frequency of index sequences of fastq file

I actually went back to the output folder on the MiniSeq and found this file: DemultiplexSummaryF1L1, which contains demultiplexing information of the whole run. Half-way through this file I found what I was looking for:
### Most Popular Index Sequences

### Columns: Sequence ReverseComplement HitCount

Index
TGACCA  TGGTCA  2455964
CGATGT  ACATCG  2320579
GGCCCC  GGGGCC  2023652
CGCGGG  CCCGCG  1983605
CCGTCC  GGACGG  1641274
CAGATC  GATCTG  1514252
CCGGCC  GGCCGG  1319383
CCGCGC  GCGCGG  1109940
ACAGTG  CACTGT  891344
CCCGGG  CCCGGG  706083
TGCCCA  TGGGCA  168294
CAGCTC  GAGCTG  164091
.... [data continues]


Tracing and quantifying inheritance over many generations

Since it seems my comment helped, I've been suggested to post this as an answer...
I'd reccomend everyone interested in these sort of population genetics to have a look at 'The Ancestors Tale' by Richard Dawkins. In particular, chapter 0, which deals with human/hominid genetics. That should help answer questions. The relevant part of the book appears to be available online for free in GoogleBooks (https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ).
Below I'll note down the answers from the comments section, not because I wrote them but I figured it might help anyone who stumbled over this question afterwards:

The chances of being a direct descendant of 1 individual of the original 10,000 is 100%, assuming they have any descendants. - Tom Melly

If one had a DNA sample of that individual, one could use uniparental markers (either Y chromosome for men, or mtDNA for women) to infer whether they are a direct descendant. - 4galaxy7

(My guess is that it doesn't have to be uniparental markers. Other genes might work just as well. I'm thinking of Neanderthal/Denisovan work. But I might be wrong.)


Hisat2 : which option should mention for strand specific library read

--fr/--rf/--ff should rarely be set, since they refer to the relative orientation of reads and basically everything these days is --fr. Mate-pairs were --rf, but you don't see those much any more.

--rna-strandedness is a very different option, since it sets how reads are expected to align against genes. Most stranded protocols in use these days follow the dUTP-method, where read #2 in a pair has the same orientation as the transcript from which it arose. So either R or RF would typically be appropriate, unless the library is unstranded. In practice, I expect this is more useful if you plan to run stringTie downstream, since then the XS auxiliary tag is set appropriately.


OpenAPI-Generator Python Client Instructions Examples

OK fine, found useful information on READme.md generated by OpenAPI-generator itself
and tried pdoc to get html document of my OpenAPI-generator generated client Module
Moreover on stackoverflow How to generate API Documentation using openapi-generator? you can find how OpenApi Generator can do generation of html documentation from your yaml or json definition of your API.
For openAPI json api specification I found RapiPDF too : Web Component to generate PDF from OpenAPI Spec (json yml file).


Best data repository to publish a large 'plotly' table containing all annotations on a transcriptome?

You may try Code Ocean. Currently, it is not that easy to use IMHO, but one of their tech leads promises to substantially improve it towards the year end.


How to deal with zeros in the denominator when calculating percent enrichment?

Cube root transformation is another way to deal with zeros without dropping them.

var^(1/3)

EDIT as Michael pointed out, square root transformation is another--perhaps more common--transformation that leaves zeros as zeros.

sqrt(var)


M <- SetIdent(M, value = "status")
or more explicitly
M <- SetIdent(M, value = M@meta.data$status)
You can also use the group.by argument of UMAPPlot() or other plotting functions from Seurat for that matter.


Is it possible to predict protein-ligand binding kinetics using machine learning?

This is a comment to redirect to other answers, but go too long.

Docking is the in silico prediction of where a ligand binds in a protein.
There are several Q&As about docking here that are relevant, so please check out:

machine learning, docking and the factors that affect it: Generate ligands candidates based on protein shape —also worth a read the "learning from the ligand" paper paper, which discusses why some ligand only parameters tell a lot about binding.
correlation between docking scores and enzyme kinetics: Purpose of correlating docking score with catalytic efficiency(Kcat/km) —discusses a paper https://www.nature.com/articles/ncomms12965 wherein a linear regression is done between docking scores and kcat/KM


Where is an up to date miRNA database and what happened to miRBase?

If you are still interested, last year miRBase generated new updates. Currently, according to ftp site the last release is 22.1. So, it is not a dead project and for more specific information you should reference the miRBase blog.
In order to get the best set of miRNA annotations, you have to define your target species or even if exists some related annotated species and then, retrieve the current annotation along different miRNA databases, but most importantly you have to be really clear about the definition of miRNA families, because it is clearly different in miRBase and also in RFAM.


How do I create a VCF file of all known pathogenic mutations in a gene of interest?

I'm not sure how to generate the additional mutations, but I would say that HGMD is not the way to find all the pathogenic variants. I would probably filter this table by either Clinical significance or Evidence->Phenotype association.


how to interpret of "pvclust" dendrogram and finding height for cutting dendrogram?

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

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

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

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

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

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


Possibility to save output blastn table in memory using biopython

You can specify send to stdout using out='-' in the Biopython wrapper.

from Bio.Blast.Applications import NcbiblastnCommandline
import pandas as pd

cline = NcbiblastnCommandline(query='seq.fna', subject='seq2.fna', outfmt=6, out='-')
output = cline()[0].strip()

rows = [line.split() for line in output.splitlines()]
cols = ['qseqid', 'sseqid', 'pident', 'length',
        'mismatch', 'gapopen', 'qstart', 'qend',
        'sstart', 'send', 'evalue', 'bitscore']

data_types = {'pident': float, 'length': int, 'mismatch': int,
              'gapopen': int, 'qstart': int, 'qend': int,
              'sstart': int, 'send': int, 'evalue': float,
              'bitscore': float}

df = pd.DataFrame(rows, columns=cols).astype(data_types)


Recommended coverage depth for ChIP Seq (H3K9me3 )

In the initial description for "usable fragments", ENCODE has the following (emphasis mine):

Usable fragment – A fragment is defined as the sequencing output corresponding to one location in the genome. If single-ended sequencing is performed, one read is considered a fragment. If paired-ended sequencing is performed, one pair of reads is considered a fragment. Fragments are considered usable if they pass the various filters in the ChIP-seq uniform processing pipelines. Used to evaluate ChIP-seq data.

Aligning reads to repetitive regions is difficult because they don't give unique alignments.
Reads from these regions will often have poor alignment scores and may not pass alignment QC checks.
Reads may also be counted as multi-mapped, and often methods call for uniquely mapped reads, which again, may not be true for a large number of reads from H3K9me3 ChIP-seq.
So for this particular histone mark, ENCODE loosens the definition for what you should count towards a minimum sequencing depth threshold.


Convert bam file to highly compressible bam

You can just set the fields you don't need to *:

samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam -

This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields.

The resulting BAM file looks like:

*   83  chr1    3003308 37  51M =   3003157 -202    *   *
*   99  chr1    3003364 33  51M =   3003547 234 *   *
*   99  chr1    3003410 37  50M =   3003636 269 *   *
*   99  chr1    3003438 255 50M =   3003760 373 *   *


AND or OR gene - protein logical relation in KEGG

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

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


This looks like genotyping failed for some samples. It's unclear from the op whether the genotyping data was controlled for quality, but I would recommend performing sample QC.
On a side note: in human gwas with imputed snps, I would expect to see clusters of significant snps as a sign that the significant snps are not an artifact. In the op there are some single significant snps. This might be a sign that those significant snps are artifacts. There are also some snps that form clusters and therefore look more promising.


What are the differences between GWAS and QTL mapping?

The basic difference between GWAS and QTL mapping is that GWAS studies the association between alleles and and a binary trait, such as being a sufferer of a disease, while QTL analysis deals with the contribution of a locus to variation in continuous trait like height.

Thus GWAS normally deals with frequencies (i.e. is the frequency disease probands greater in those with the minor allele compared to the major allele at locus X), whereas QTLs deal in correlations (Is the number of minor alleles at locus Y correlated with height).

Methodologically the primary difference would be the distribution that links the allele to the outcome (disese/continuous phenotype). GWAS is going to use Chi-squared statistics or logistic regression, QTL with more likely use straight up normal distribution based linear regression.


combining the use of workdir and option --jobscript in snakemake

To quote from the snakemake documentation

All paths in the snakefile are interpreted relative to the directory snakemake is executed in. This behaviour can be overridden by specifying a workdir in the snakefile:
workdir: "path/to/workdir"
Usually, it is preferred to only set the working directory via the command line, because above directive limits the portability of Snakemake workflows.

It's not explicitly stated, but it's sort of implied that all relative paths then become relative to the working directory. I would expect that specifying an absolute path would get around this.
As an aside, in my mind setting the working directory is usually only needed on clusters without a shared file system (presumably shared between the worker nodes, but not with the head node), since there you can't cd $WORKDIR before running snakemake. This is then normally done with your scheduler, in such cases.


Overlap in Paired-end Reads for Sequencing?

Yes, if the insert size is smaller than the read sizes, this would happen. For some applications, for example SNP detection for molecular diagnostics purposes, this approach is used.


Low Fraction of usable antibody reads in CiteSeq

We have also seen a very similar issue and are also investigating why. Possibly of note, I tried forcing the cellranger program to assign a very large number of cells (10x the actual number of cells we had) using the --force cells option. This appeared to increase the % of usable Ab reads and Ab reads per cell metrics, suggesting that these reads may be associated with barcodes that have relatively few UMIs?


Best distance parameter for estimating physical interaction between residues in a PDB file

A popular (and I would say, respected and trusted) website is PDBsum https://www.ebi.ac.uk/pdbsum (which also has a Wikipedia article about it: https://en.wikipedia.org/wiki/PDBsum)

They measure protein-protein (prot-prot) contacts as any Nitrogen, Carbon, or Oxygen element atoms (N, C, O) within 4.0 Angstroms (Å) when measuring the 3d euclidean distance.

If your atomic contacts are not close enough, the amino acids may not actually interact at all, even though they theoretically could.  Since in PPIs you generally want to use reliable ground-truth information, it is better to use cautious parameters.  Often, in PPIs, the answer to "does it interact?" is not "yes" or "no" but rather "yes" or "don't know / unknown".


Why not use sed instead?

sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

Or even simpler:

sed 's/chr_//' mySequence.fasta > mySeq.fasta


The visualisation of a list of genes on URD object

I just added URD clustering labels for each cel to seurat object and then I did DoHeatmap by using group.by="URD clustering labels"


Searching for gene expression data by cell line

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


How did researchers derive the Ramachandran "validation" contours?

To expand on Matteo Ferla's comment as an answer, there are two approaches to getting these "allowed" contours:

From probabilities found in observed crystal structures. You could imagine plotting a Ramachandran plot for the whole PDB and colouring the points by assigned secondary structure from DSSP. The contours would then be regions that encompass e.g. 99% of alpha-helical residues, or 99% of all residues. In doing this you would notice for example that glycine can adopt more regions of the Ramachandran space that other amino acids. PROCHECK is software that scores a protein based partly on the likelihood of the phi/psi angles of each residue.
From reasoning about the chemistry of the secondary structure of proteins. Only certain phi/psi combinations allow the i -> i + 4 pattern in an alpha-helix or the alternating up/down sidechain pattern in a beta-strand.


Extracting features or gene from PCA after calculating PCA for downstream analysis

The help page explains that there is a component which:

  rotation
  the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:

C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.


You make a phylogenetic tree of the "host species" and compare it to the phylogeny of the TE's. Homology is not a good idea for evolutionary trees.

There are two snags:

Concerted evolution, the TE's will risk recombining with each other within a host genome and that can mess up the phylogeneic signal (the theory as why is quite complicated)
If host species phylogeny are higher eukaryotes thats a lot of genetic material, so you'd need a short-cut. You can either use someones species tree, or focus on TE's within a subsection of the host genome and produce a tree of this subsection.

With point number 2 "a subsection of the genome" you can then assess synteny, i.e. has the TE jumped between position between species. If the TE hasn't moved then its an orthogonal proof of putative vertical transmission.

What you are looking for is incongruence between host phylogeny and TE phylogeny as evidence of horizontal transmission and congruence to call vertical transmission. However, incongruence could occur due to the complexity of concerted evolution, so synteny would assist again as an orthogonal proof. Congruence is assessed as a "mirror tree", there is quite a lot of theory on this looking at host-parasite coevolution

For making trees RAxML or PhyML are flavours of the month particularly for TEs.

A state of the art-of-the-art review is here, but it goes into significant detail into precise proofs, which in my opinion arn't needed


What is the best way to start with Structural bioinformatics with zero prior experience?

As others have said, there is no short cut and you will need to spend some time with a good book.

I would say there are 3 broad areas of understanding you have to have as a structural bioinformatician:

Solid understanding of the basic science. This includes molecular biology but also chemistry, which is often overlooked.
Knowledge of the databases and data formats available in the field. You need to know where to go to get protein structures, protein sequences etc.
Programming and algorithmic skills. You don't have to be a computer whiz but being effective in a high-level language (e.g. Python, R, Julia) and being able to parse data and carry out basic data analyses are essential.

If you want to take it further and have a bit more time then I would consider implementing something yourself such as a protein secondary structure predictor. Getting this working is challenging and would show you have competency in the above three areas.


Proteomics data Vs Transcriptomics data?

This is what I found from doing some research. Comments are welcome at any point!

Capture percentage for data gathering:

Transcriptomics Data: There's an amplification step in Transcriptomics data gathering methods. Hence, it's possible to capture almost the totality of the Transcriptome using those methods(scRNA seq methods, Nanopore tech, Spatial transcriptomics, etc)
Proteomics Data: There's no amplification step in Proteomics data gathering methods. Those methods have fragmentation and a rejoining step which leads to a loss of the data.

Scalability of methods:

Transcriptomics Data: Scalable methods available
Proteomics Data: Less scalable methods for protein studies

Reference availability:

Transcriptomics Data: Fully annotated references available on consortiums such as ensembl biomart, etc
Proteomics Data: Universal and Comprehensive Human Proteome reference is still in question and because of the incomplete and inaccurate references, much of the newly generated data is being rendered useless

Uniformity

Transcriptomics Data: Lots of publications on pipelines, and their benchmarking available
Proteomics Data: Lack of uniformity across labs/research groups and lack of related literature lead to differences in protein fragmentation and solubilization and differences in algorithms to run analyses

Technical Bias

Transcriptomics Data: Many methods available to tackle cell-bias, noise in data, and batch effects, so as to capture maximum biological variability
Proteomics Data: Cell-wide Mass spectrometry has a bias towards the identification of peptides with higher concentration or contamination from other experiments


what percentage of the human genome is MAPQ=0?

Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper is not necessarily MAPQ=0 to another.

That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a useful answer. For 150bp reads, each reference position is covered by 150 reads. What if 50 of them have no other exact hits elsewhere, but the rest 100 have? Is this position a repeat or not? In addition, what if the 50 are unique only because one mismatch? If there is a variant at that mismatch, the 50 would become repeats or mapped elsewhere.

My preference is to say a position is "unique" under k-long reads if over k/2 reads overlapping the position have no other perfect or 1-mismatch/1-gap hits elsewhere in the genome. Under this definition, 79.3% of human genome are unique for 35bp reads. 92.4% for 75bp reads. I don't have the number for 150bp reads. I guess it will be around 95%. Empirically, 94-95% of human genome is callable with 100bp paired-end reads.

As to other measurements, the most common one is the fraction of reads that has an exact hit elsewhere. You can use Fabio's method. It gives a good enough estimate once you simulate over 1 million reads. This fraction is around 86% for 35bp reads and 95% for 75bp reads, as I remember. The problem with this approach is this fraction is not very informative to variant calling due to the issues I talked about. Another way is to use RepeatMasker. It is worse. RepeatMasker masks 50% of human genome, but excludes segmental duplications where short reads can't be confidently mapped.


The snarky answer is that it's better to find a matching dataset at the experimental design phase (i.e. before these studies are carried out), so that the bioinformatics component on the other side of things ends up being a lot easier and more robust. Although it can be quite time consuming to do a database search at the start, it's even harder to find matching projects if the legwork hasn't been done during the design phase. Researchers should be thinking about how hypotheses will be validated before the experiments are carried out.

Every additional characteristic that is chosen for comparison adds complications with regards to matching data, and reduces the applicability of results when there is a mismatch. You would have to be somewhat lucky to be able to find a IBD-Crohn's dataset with paired RNA-seq data and 16s-seq (presumably short-read; is 454 good enough?) from a specific non-stool tissue.

Apart from the places you've already looked, there's dbGaP and SRA. Given the number of characteristics that have been mentioned, I expect you'll have to make some compromises when comparing other datasets.


processing multiple fastq files with cutadapt

Try this one

for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG "$file" >> results.out
done

Also you have omitted --discard-untrimmed from the original command. If you want that as well along with everything to run in parallel:

for file in /dir/*
do
  cutadapt -g ACTTAAGTGTATGTAAACTTCCGACTTCAACTG --discard-untrimmed "$file" -o "`basename -s .fastq $file`_out.fastq" &
done


Identify differentially covered genes only between two samples

Not really an answer but an extended comment... and most likely something you don't like to hear

I guess by technical replicates, it means taking the same biological sample and making 2 methylation libraries. If this is used as replicates in deseq2 or edger, the variance you are estimating is the technical variation that comes with preparing the library, bisulphite conversion etc..

Often the biological variance will be much larger. That is if you take 2 independent samples of group A, their variation in methylation is much larger than the technical variation you see from handling the samples. If the technical replicates are treated as biological replicates and differential analysis is performed, it will be seriously misleading.

I think both edgeR and DESeq2 doesn't estimate dispersion with 1 v 1 (which is right).

I would suggest maybe the following:

Run a fisher.test at every gene (methylated / unmethylated), get the odds ratio and p-value, rank this, but explain that the p-value is only for ranking, and will underestimate the variation.
A bit more refined, look for other studies that might be close to yours, get an estimate of the dispersion in their study and use a negative binomial distribution to again rank genes that show the biggest difference


By default, DimPlot colors the cells by their "identity class", stored in Ident(scData) in your case. This indeed contains the clustering information.

However, you can you use the group.by argument for DimPlot to pass another column from scData@meta.data to color cells by. In your case you would need a metadata column that stores the ptx/human information.

The merge function automatically stores the original "project" information of the individual Seurat objects that you combine under orig.ident in the metadata. Since you set the project names to either "ptx" or "human", the orig.ident column should contain the information you need.

Something like this should work:

DimPlot(scData, reduction = "umap", group.by = "orig.ident", label = TRUE)


Get start and end coordinates per chromosome

Not sure if I completely understood your question, but if you want simply an interval of your entire genome  you could do instead:

samtools faidx MyGenome.fasta

Obviously that means that you need access to that, not sure if you have.
The resulting fai file contains a format similar to that:

000000F 33203223        94      60      61
000001F 28828106        33756799        60      61
000002F 27810542        63065468        60      61

Where the first column is your chromosome and the second one your length.
I guess you could then just do a awk '{OFS=FS="\t"}{print $1,"0",$2-1}' to get your file.


As @terdon mentioned in the comments, you need to sudo yum install libicu on the relevant nodes (or have your admin do that for you). Alternatively, seurat can be installed via conda, which means you don't need root access.