Differential gene expression analysis of time series with replicates

R packages like DESeq2, edgeR, and sleuth will help you perform differential gene expression analysis.
These are typically used for control-vs-treatment experimental designs.
So you could use them like this and perform a bunch of pair-wise comparisons between the time points (i.e. what are the differentially expressed genes between timepoint $t_0$ and timepoint $t_1$, $t_0$ and $t_2$, etc.).
But all of these methods make use of generalized linear models to perform the statistical inference on coefficients.
For the simple control-vs-treatment design, above, this looks something like
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i X$$
where $g$ is the link function, $\mu_i$ is the mean expression of gene $i$, $\beta_i$ is the expression fold change estimate, and $X$ is the indicator variable saying if the sample is a treatment or control sample.
You can modify this equation by adding a time component to it.
$$g \left(\mathbb{E} \left[Y_i | X \right] \right) = \mu_i + \beta_i^{(group)} X + \beta_i^{(time)} t$$
If you specify the timepoint of the sample in the experimental design matrix for DESeq2, edgeR, or Sleuth, you can get information about the effect of the perturbation as well as the effect over time.
You can try more complicated designs, too, like a coefficient for the interaction between timepoint and perturbation status, but with 3 replicates in each, it's probably best to stick with simple, linear models.


Are codons in RNA layered? Are we misinterpreting RNA codons?

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


If you change the delimiter to something besides a space in the tr command (tr "\"" "\t") then you can use that same delimiter in the awk command (-F "\t").  In my example below I used a tab, also I filtered the taxon prior to the tr so that you don't have to worry about that staying in sync if you change the delimiter.

 egrep -v "/db_xref=\"taxon"| tr "\"" "\t" | egrep "/gene=|/db_xref=|/product=|/protein_id=" | sort | uniq | awk -F "\t" '{print($2)}END{print("")}';

results:

GeneID:100933602
UROS
uroporphyrinogen-III synthase isoform X1
XP_003755055.2


Source for whole genome comparisons at NCBI Genomes

I inquired into the details of the dendrograms after becoming frustrated with the lack of information. As with Ensembl, I'm sure that the folks at NCBI have a standardized pipeline that they run the sequences through to generate these dendrograms as no specific source deserving attribution was mentioned:

The tree is based on a pairwise, BLAST comparison of chromosome sequences (an assembly is used for genomes with 2 or more chromosomes). BLAST identity is normalized on sequence length. The distance map is generated by a Neighbor Joining algorithm, where the distance is the blast score.
Dr. Wayne Matten, NCBI


scater: SingleCellExperiment, Error in seq_len(ncol(assay))

I think you might want to use the most recent workflow, which should be for the latest version on Bioconductor.


Intersecting two different files with one "master" file based on different columns


Store contigs and genome substrings based on mismatches

For the first question try something like:

new_substrings = []
new_positions = []
new_list = []
for pos, sub in zip(positions, substrings):
   if not pos:
       new_list.append(sub)
   else:
       new_substrings.append(sub)
       new_positions.append(pos)

this gives:

>>> new_substrings
['TACA']
>>> new_positions
['11, 4']
>>> new_list
['CTAGAAAAG']

For the hamming distance you just need to iterate over the values in seq_ref and append the results of the hamming distance function in mismatches.  Based on the given code and the example above you should be able to create that block.  Below is an untested quick example to get you on track.

mismatches = []
for seq in seq_ref.keys():
   for value in seq_ref[seq]:
       mismatches.append(HammingDistance(seq, value))

As far as the CSV and sorting is concerned I suggest that you look into pandas or numpy.  Both of those packages provided functions to help with this.


Try LR_Gapcloser.

I've used L_RNA_Scaffolder for trying to scaffold a genome (which turned out to be more complex than I had expected). It looks like LR_Gapcloser (written by the same people) is similar, but designed specifically for scaffolding using long-reads.

That page also suggests PBJelly and GMCloser as competing tools.


What are the ways to process a list of differentially expressed genes?

You originally had asked a very broad question, so I'll try to demonstrate why that is such a hard question to answer. I've done two fairly large differential analysis studies (and a few smaller ones) covering very different areas of research, and the approaches that other researchers used subsequent to my differential expression calculations were unsurprisingly also very different.

In the first study, I did a de-novo transcriptome assembly of Schmidtea mediterranea (Smed) from Illumina RNASeq reads, then carried out a DE analysis (DESeq) on a β-catenin knockout mutant population compared to wildtype. The other researchers had some very impressive tools available at their disposal, allowing them to subsequently visualise the positional expression of differentially-expressed genes in the knockout and wildtype populations after amputation. They found a few candidate genes that had graded expression along the anterior-posterior axis in the wildtype population, and this expression disappeared in the knockout mutants. One of these genes was a teashirt (tsh) family member, and a good candidate for further discovery based on the researchers' existing knowledge of the gene pathways associated with cellular regeneration. To cut a very long story short, they experimentally validated the teashirt link through further RNAi studies, generating a double-headed Smed as a result. This was followed up in Zebrafish (Danio rerio), because there happened to be another building about 100m away that specialised in Zebrafish development.

The second study was/is a bit more similar to your current situation. We looked at sorted dendritic cell expression in mice after the application of two Th2-inducing treatment conditions (injected Nippostrongylus brasiliensis; applied DBP+FITC) and one control condition (injected PBS). I did a differential expression analysis with DESeq2, and developed a Shiny App for result exploration after I got tired of getting asked to do a whole bunch of mechanical follow-up tasks (mostly generating graphs, lists, and heatmaps for different gene subsets). We ended up with thousands of differentially expressed genes, so it wasn't feasible to look individually at every gene for research evidence and pathways -- one scientist spent a couple of weeks looking at one interesting-looking candidate gene. The differentially-expressed genes were fed into the oh-so-expensive IPA (which gives results, but not necessarily good ones) and DAVID, among a few others. These helped a little bit, and led to what was the central idea of our paper (i.e. that there were at least two different Th2 responses), but we still had the problem of too many "significant" pathways, and lots of differentially-expressed genes that were members of a whole bunch of other pathways. We tried WGCNA; we did a bit of GSEA; and we compared Limma and DESeq2, but still couldn't find anything concrete or obvious to demonstrate what was triggering the Th2 response. There were experimental studies as well: blocking and mutant studies to investigate the two types of Th2 response (and experimentally confirm the differential expression). Our experiments and exploration are still ongoing, but we felt it was necessary to get one paper out at the intermediate stage, so just ended up taking the top 25 differentially-expressed genes in each experiment (by fold change), and finding some way of grouping some of them together (via expert knowledge and literature searches). It kind-of-sort-of worked, but we're left with a bit of frustration in that we (and the tools we used) were basically stumped by having too many biologically-relevant results.


how can i download rna-seq data sample for starting rna-seq analysis by R?

You must have found the airway dataset by googling "RNASeq R tutorial"

https://www.bioconductor.org/packages/devel/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html

Why won't this work for you?


Error using htseq-count: Could not retrieve index file

The problem seems to be missing index files of the bam files (therefore the error message that the index file is missing). You should be able to generate those with samtools, e.g.
samtools index ./OVCA429_587680_1_1.rnaseq.fastqAligned.sortedByCoord.out.bam


How to subset genes and its nested features from a GFF file using a gene list

Given the nested structure of a GFF file, it may be easier to make a python script using a GFF file parser, like gff3, gffutils, or BioPython.
Dealing with nested data in bash or with some fancy awk script is probably going to trickier to write than it's worth.
You can start by ensuring your list of target genes is sorted in the same order as the GFF file, then iterate through the file. When you encounter a like with mRNA in the 3rd column, store that ID value and extract all the subsequent lines with that same ID.
If the line does not have that ID, you know you've collected all the genes, mRNAs, etc associated with that gene and can pop it from your queue.
Repeat until you reach the end of the GFF file or your list of target genes is empty.


You can use a list and insert() to add an element in a specific order, then expand the list with *. Or you can use join().
from Bio import SeqIO
import sys

for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):
  sequences = []
  sequences.insert(0, hello_fasta.description)
  sequences.insert(1, len(hello_fasta.seq))
  sequences.insert(2, hello_fasta.seq)
  # option 1
  print(*sequences, sep='\t')
  # option 2
  print('\t'.join(map(str, sequences)))


Map gene IDs to protein structure if exist

Not sure if you want the structure exactly or a handle of some sort, but the best way i know to get ALL your 5000 identifiers to map to uniprot/pdb is uniprot's mapping service.
You can either copy your ids straight from excel or parse with python ex.:

import pandas as pd
import os

file_name = 'aap7847_Zhang_SM_Table-S5.xlsx'
sheetname = 'Table S5'
df = pd.read_excel(file_name, sheet_name=sheetname, skiprows=1)
ids = list(df['Gene_ID'])

name_of_file_to_save_to = "geneIds.txt"
with open(name_of_file_to_save_to, "w") as newf:
    for id in ids:
        print("{}\n".format(id), file=newf)
    newf.close()

You can then save the Uniprot's response and either add it as a column to your file or map to other parts of Elixir ecosystem with SIFTS, for example, as @Matteo Ferla says.

In the case that you maybe need a single id's information, you could hit their API too, although this wouldn't work for bulk or at least didn't work for me. You can map the uniprot accession you get from this to PDB or whatever else you need!

import requests
def get_uniprot_accession(geneid, format='list'):
    resp = requests.get(
        "https://www.uniprot.org/uniprot/?query={}&sort=score&format={}".format(geneid, format)).text
    if format == 'list':
        resp = resp.strip()
    return resp


Efficiently iterate over bed and subset pyranges objects

Setup:

import pandas as pd
import pyranges as pr
from pyranges import PyRanges
from scipy.stats import fisher_exact
import numpy as np

# ! zcat dataset1.tsv.gz | head -2
# chromosome    start   end num_motifs_in_group called_sites    called_sites_methylated methylated_frequency    group_sequence
# chr21 5010053 5010053 1   3   0   0.000   CACCACGTCCA

# avoid using more columns than we need, especially seq which contains object
# is heavy to lug around
colnames = ["Chromosome", "Start", "End", "calls", "methylated"]
usecols = [0, 1, 2, 4, 5]

h1 = PyRanges(pd.read_csv("dataset1.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))
h2 = PyRanges(pd.read_csv("dataset2.tsv.gz", sep="\t", names=colnames, header=0, usecols=usecols))

regions = pr.read_bed("chr21.bed.gz")
regions = regions[~regions.Chromosome.str.contains("_")].merge()

regions.ID = list(range(len(regions)))

regions
# +--------------+-----------+-----------+-----------+
# | Chromosome   | Start     | End       | ID        |
# | (category)   | (int32)   | (int32)   | (int64)   |
# |--------------+-----------+-----------+-----------|
# | chr21        | 5065660   | 5065810   | 0         |
# | chr21        | 5068140   | 5068290   | 1         |
# | chr21        | 5128400   | 5128550   | 2         |
# | chr21        | 5154780   | 5154930   | 3         |
# | ...          | ...       | ...       | ...       |
# | chr21        | 46637040  | 46637190  | 2377      |
# | chr21        | 46638480  | 46638630  | 2378      |
# | chr21        | 46639820  | 46639970  | 2379      |
# | chr21        | 46692680  | 46692830  | 2380      |
# +--------------+-----------+-----------+-----------+
# Unstranded PyRanges object has 2,381 rows and 4 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

h1
# +--------------+-----------+-----------+-----------+--------------+
# | Chromosome   | Start     | End       | calls     | methylated   |
# | (category)   | (int32)   | (int32)   | (int64)   | (int64)      |
# |--------------+-----------+-----------+-----------+--------------|
# | chr21        | 5010053   | 5010053   | 3         | 0            |
# | chr21        | 5010215   | 5010215   | 5         | 1            |
# | chr21        | 5010331   | 5010335   | 6         | 0            |
# | chr21        | 5010525   | 5010525   | 2         | 2            |
# | ...          | ...       | ...       | ...       | ...          |
# | chr21        | 46699713  | 46699718  | 274       | 254          |
# | chr21        | 46699739  | 46699739  | 133       | 113          |
# | chr21        | 46699756  | 46699756  | 37        | 32           |
# | chr21        | 46699770  | 46699773  | 156       | 130          |
# +--------------+-----------+-----------+-----------+--------------+
# Unstranded PyRanges object has 331,373 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Get total number of methylated per region:

def count_total_and_methylated(df):
    """For each region in bed, get sum of calls and methylated"""

    grpby = df.groupby("ID")
    total = grpby.calls.sum()
    methylated = grpby.methylated.sum()

    chromosome = df.Chromosome.iloc[0]
    start = grpby.Start.first()
    end = grpby.End.first()
    region_id = grpby.ID.first()

    return pd.DataFrame({"Chromosome": chromosome, "Start": start, "End": end, "ID": region_id,
                         "calls": total, "methylated": methylated})

rh1 = regions.join(h1).drop(like="_b").apply(count_total_and_methylated)
rh2 = regions.join(h2).drop(like="_b").apply(count_total_and_methylated)

rh1.File = 1
rh2.File = 2

def add_missing(df, df2):

    """Since join is inner rh1 and rh2 might not have the exact same intervals. Fix this."""

    missing_in_first = np.setdiff1d(df2.ID, df.ID)
    first_missing = df2[df2.ID.isin(missing_in_first)].copy()

    f = df.File.head().iloc[0]

    first_missing.loc[:, "File"] = f
    first_missing.loc[:, "calls"] = 0
    first_missing.loc[:, "methylated"] = 0

    return pd.concat([df, first_missing])

rh1 = rh1.apply_pair(rh2, add_missing)
rh2 = rh2.apply_pair(rh1, add_missing)

j = rh1.join(rh2, suffix=2).drop(like="File*|ID*|(Start|End)2")
j
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# | Chromosome   | Start     | End       | calls     | methylated   | calls2    | methylated2   |
# | (object)     | (int32)   | (int32)   | (int64)   | (int64)      | (int64)   | (int64)       |
# |--------------+-----------+-----------+-----------+--------------+-----------+---------------|
# | chr21        | 5065660   | 5065810   | 188       | 8            | 70        | 0             |
# | chr21        | 5068140   | 5068290   | 168       | 166          | 74        | 72            |
# | chr21        | 5128400   | 5128550   | 901       | 14           | 143       | 1             |
# | chr21        | 5154780   | 5154930   | 653       | 0            | 593       | 0             |
# | ...          | ...       | ...       | ...       | ...          | ...       | ...           |
# | chr21        | 46637040  | 46637190  | 66        | 63           | 69        | 62            |
# | chr21        | 46638480  | 46638630  | 28        | 0            | 30        | 3             |
# | chr21        | 46639820  | 46639970  | 22        | 0            | 18        | 14            |
# | chr21        | 46692680  | 46692830  | 71        | 25           | 57        | 23            |
# +--------------+-----------+-----------+-----------+--------------+-----------+---------------+
# Unstranded PyRanges object has 1,922 rows and 7 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Compute fisher exact using fisher instead of scipy:

def fe(df):
    from fisher import pvalue_npy

    cols = "calls methylated calls2 methylated2".split()
    df[cols] = df[cols].astype(np.uint)
    _, _, two_sided = pvalue_npy(df.calls.values, df.methylated.values, df.calls2.values, df.methylated2.values)

    t1, t2, m1, m2 = df.calls, df.calls2, df.methylated, df.methylated2

    _or =   (t1 * m2) / (t2 * m1)

    df.insert(df.shape[1], "PValue", two_sided)
    df.insert(df.shape[1], "OR", _or)

    return df

result = j.apply(fe, nb_cpu=1)
result.drop(like="methylated*|calls*")
# +--------------+-----------+-----------+------------------------+--------------------+
# | Chromosome   | Start     | End       | PValue                 | OR                 |
# | (object)     | (int32)   | (int32)   | (float64)              | (float64)          |
# |--------------+-----------+-----------+------------------------+--------------------|
# | chr21        | 5065660   | 5065810   | 0.1152954430992567     | 0.0                |
# | chr21        | 5068140   | 5068290   | 0.9999999999995065     | 0.9846955389124064 |
# | chr21        | 5128400   | 5128550   | 0.7075919513120225     | 0.4500499500499501 |
# | chr21        | 5154780   | 5154930   | 1.0                    | nan                |
# | ...          | ...       | ...       | ...                    | ...                |
# | chr21        | 46638480  | 46638630  | 0.24262295081967644    | inf                |
# | chr21        | 46639820  | 46639970  | 0.00024909594620266284 | inf                |
# | chr21        | 46692680  | 46692830  | 0.7354565511951735     | 1.1459649122807019 |
# | chr21        | 36784100  | 36784250  | 1.0                    | nan                |
# +--------------+-----------+-----------+------------------------+--------------------+
# Unstranded PyRanges object has 1,932 rows and 5 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.


creating graph of distribution of blast hits on the query sequence

You probably want to include query start (qstart) and query end (qend) in your blast output.

Something like this:

blastn -query your.fasta -out blast.out.txt -db your.db -outfmt '6 qseqid sseqid qstart qend length evalue'

In R you can take the "qstart:qend" from each line for density plot.

There are many ways in R to plot the densities of these start and end amino acids.

Let me show an example with a small data frame:

qstart <- c(200, 300, 250, 400, 500)
qend <- c(300, 450, 400, 600, 650)

df <- as.data.frame(cbind(qstart,qend))

aa <- vector()
i=1
for(i in 1:5){
aa <- append(aa, c(df[i,1]:df[i,2]))
i+1
}

hist(aa)
dens <- density(aa)
plot(dens)


Should PCA be standardized for gene expression?

When the gene expression is scaled and centred you reduce the difference between genes.

Imagine you have one gene A that is highly expressed usually and has a standard deviation of 500 units compared to a gene B that is not much expressed and only have a standard deviation of 5.

In the scaled and centred genes both contribute the same because A usually is expressed 10000 and B is usually expressed 100 units. So, for both the standard deviation is 5% of their expression. Meaning that a variation in one is as important as in the other.

If not scaled (and centred) the first gene A will contribute more to the variation than gene B, because the expression variation in absolute numbers is bigger.

Both are used in publications, but I think the scaled and centred is more used, because the first dimension reflects "better" the differences between samples. Of course you can do both and select the one that is better to show what you want.


Does SBOL support timing and threshold value parameters?

Note: I will use SBOL3 in my descriptions; the answer can be equivalently translated to SBOL2
At present, SBOL doesn't have any explicit recommendations for how to represent these types of timing parameters. This is because there is not yet a clear scientific consensus on whether or not such parameters are a good way to represent genetic circuits.
That said, one of the advantages of using SBOL is that its extensible RDF-based model allows additional information like this to be added into a model. Two good ways of doing so are:

If the parameters are simple and can stand alone as measurements, then they can be added to designs as a Measure object with appropriate units and type. For example, if you have biological devices with a 90 minute propagation delay, then on the Component representing the device you can add a Measure with numerical value 90, units of om:minute.  You'd likely want two types, one of SBO:0000225 (Delay), and a custom type for your own notion of propagation delay, since I don't think that's included in any ontology yet.
If the parameters are more complex and/or inter-related, e.g., a full data-sheet, then you'd likely want to instead create one or more extension classes that are used to represent the information. A library that makes this easy is sbol_factory, though one can do it with any of the SBOL libraries.

This becomes more complicated if you want to talk about behavior in different environmental contexts, but the basic principles are still the same.
Finally, as for SBOL visual: again, there's no explicit representation, but the model is extensible, and this information can be added as any sort of text or custom glyph as long as it doesn't conflict with ones already in the standard.


Deciding which samples go in which batch

I found the OSAT package which does almost what I want.

However it is specific for Illumina and for a limited number of plates: "plates that hold 2, 4, 8 IlluminaBeadChip chips and have 24, 48, 96 wells, respectively."
It doesn't take in consideration the numeric values:

library(OSAT)
inPath <- system.file('extdata', package='OSAT')
pheno <- read.table(file.path(inPath, 'samples.txt'),
                    header=TRUE, sep="\t", colClasses="factor")
pheno$AgeGrp <- runif(576, 0, 100)
gs <- setup.sample(pheno, optimal=c("SampleType", "Race", "AgeGrp"))
gc <- setup.container(IlluminaBeadChip96Plate, 6, batch='plates')
gSetup <- create.optimized.setup(sample=gs, container=gc, nSim=1000)

The minDiff package reported in Stats.SE.

It handles both numeric and categorical data, but it can only optimize for two nominal criteria.

I wrote experDesign which does everything I wanted

It can handle both numeric and categorical data for an arbitrary number of variables. It is not restricted to an even number of samples per batch.


DMPfold: number of models produced using "10 500" feature in "run_dmpfold.sh"?

The sequence_n.pdb files are temporary and are combined into an ensemble.i.pdb file after each iteration i. That file has many structures in. After the run is complete there will be a final_1.pdb file. Don't worry about the cluster files. If you are interested in what they do, consult the run_dmpfold.sh script and see what it is doing line by line.

As I have said multiple times elsewhere, you should run the script with default parameters and let it run to completion. Then use the final_1.pdb and any other final_*.pdb files as the output. If you want more structures, use the ensemble.i.pdb files.


How to subset a BAM by a list of QNAMEs?

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

  READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.


Remove/delete sequences by ID from multifasta

If you want to learn how to do things with command-line tools, you can linearize the FASTA with awk, pipe to grep to filter for items of interest named in patterns.txt, then pipe to tr to delinearize:

$ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" > out.fa

This will work on multiline FASTA.


How to check if indels in VCF files are left or right aligned?

One possibility is to run bcftools norm on the file. At the end it will print out a statistic about how many variants were realigned. I did this for dbSNPv151 (GRCh38) only for chromosome 22 like this:

$ bcftools view -i 'TYPE="indel"' dbSNP_GRCh38.p7_b151.vcf.gz 22|bcftools norm -Ou -f GRCh38.fa > /dev/null
Lines   total/split/realigned/skipped:  901503/0/490/0


Download data from the Human Microbiome Project via ascp

I can't test for ascp because they do not longer support my current setup but you can follow the FTP instructions. With FileZilla:

url: public-ftp.ihmpdcc.org


Optimum parameters for genome assembly via canu

The canu parameter reference guide is always a good place to start.

It would help to know a little bit more about what your assembly looks like. E.g.

what is your n50, n90...?
What is your dataset coverage?
Have you run any tools like quast to evaluate quality?
What is your input data type (probably PacBio)? Read length distribution?
What is your expected chromosome/replicon number?
Have you compared the assembly to any other data, e.g. genomes of close relatives? (quast can do this)

Basically, why do you think that 200 contigs is too many contigs? Sometimes you just don't have enough coverage to get more contiguous/"nice" assemblies. That's not to say it's not possible, of course, just that we don't have enough information to evaluate.

It's possible that some of those contigs are small and artefactual, so you might get rid of them with just length/coverage filtering. You can do some of this in canu itself (see "Output Filtering" section).


Is the genome a programming language (i.e. LISP)? Can we analyze it with computer science theory?

I see what the metaphor has been inspired by, where I think it got turned around, and might have a better one.

I'd like to back up a bit and establish some common ground.  It isn't wrong to expect that some aspects of genetic code and software code are similar, and while biology is usually very specific, some aspects can be seen productively from the light of computer science ideas.  For example, mass-action molecular dynamics, which are often a good model for protein interactions, are Turing-complete.  Here's a fun example where people have written a programming language out of molecular dynamics:  https://arxiv.org/abs/1809.07430

It's also true that proteins have particular substructures leading to certain functions, such as using alpha-helices to react less with the outside world.  Learning about larger substructures, called protein domains, help identify the functions of families of proteins.

The problem is that the parts of the genome which appear nested like LISP cells actually aren't included in the protein.  The regulatory regions, promoters, inhibitors, introns, and so forth never become part of the key "computational" machinery.  This nesting isn't actively part of the computational activity proteins do perform.  The main quality of LISP is that the structure and function are perfectly self-similar throughout an expression, but the structure of the regulatory regions are neither replicated down into the protein nor are necessarily recapitulated on a larger scale.

Unfortunately, we don't expect the computational structure of proteins to be the same from one protein to another.  This is because "binds only with X" and "binds only with Y" must be implemented differently to successfully be distinct operations, and it is these binding constraints by which molecular dynamics implement "computational" operations, such as "not X" and "not Y", and so those are distinct as well.

I think a better way to think about the regular structures of the regulatory apparatus is more like a network protocol.  Along an otherwise undifferentiated binary channel, one has to identify when and where the data is delivered to an application.  Similarly, along a strand of bases, the machinery of translation has to have an indication of where the "application" of the protein starts and ends, and when that protein should be delivered to the cell.  We might compare some of the images of genetic structure that provided in the original question with layouts of network headers: http://slideplayer.com/slide/4798836/15/images/16/TCP/IP+Packet+Structures.jpg  We might compare bacterial versus multi-cellular regulation as having different delivery mechanisms and data structures for delivering the application information, for different functional consequences, in the same way TCP and UDP have different structures for different trade-offs.

This is all complicated somewhat by the fact that proteins do then interact with the regulation of genes, leading to interesting "computational" dynamics such as feed-forward regulation.  However, I'm not aware of any particular structure of the regulatory apparatus that necessarily correlates with these larger functional patterns.  The opposite is frequently biologically useful: the same structure of a protein helping to activate one gene might simultaneously inhibit another.

One book I've found very accessible and useful in understanding where computation is useful biologically is Dennis Bray's "Wetware: a computer in every living cell".  It gave me a clearer picture of how computation occurs in cells.  That should open up some introductions to systems biology which explains the biological-role of larger molecular "computational" patterns.


Survival analysis using CoxPH - Effect of covariates

Example 2 and 3 seem symmetric because they are forced to by the model. The Cox proportional hazard relies on the proportional hazards assumption. This basically means that the ratio of the hazards are constant over time. If you were to plot the hazard ratio, it would look like a flat line over all Event_at. Another way of saying this is that the hazard ratios are constant over time.
How that shows up in a survival plot is a little different because of the relationship between hazards and survival. The survival curves end up following the pattern exemplified in Example 2.
As for example 3, I would guess that there are few individuals with event times past 3500. This is why you have those large jumps.
Proportional Hazards
So what does the proportional hazards assumption mean? Below are some plots of hazard functions to help your intuition on what the proportional hazards assumption means. In all of the following scenarios, the hazard ratio is 2 for blue-to-red. However, the proportional hazards assumption is only true in A, B, and C.

Since the proportion is remaining constant, it can look a little hard to see. To make it easier on the eyes, you can log-transform the hazard functions. Below is a plot that does that.

As you can see, the hazards are a constant distance from each other for A, B, and C. Basically, the proportional hazards assumption is saying that your scenario looks something like one of those plots.
In scenario D, the hazard ratio varies over time. At the beginning it is basically HR=1. At the end it is basically HR=4. With the proportional hazards assumption, Cox regression will essentially take the 'average' of the hazards over time. In this case it would give you an estimate of HR=2. Depending on the scenario, this 'smoothing' may be more or less acceptable. But this is basically what the proportional hazards assumption does in the background.
The proportional hazards assumption is actually more problematic in the cases where the two hazard functions cross each other. Below is an example of that.

In this scenario, blue has an increased hazard at the start (relative to red) then a decreased at the end. A Cox proportional hazards model would give you HR=1 in this case, because the average over all time points is no relationship.
So how do you tell if you violate the proportional hazards assumption? You may want to use a stratified Kaplan-Meier by your variables to get an idea of what the hazards look like without stipulating the proportional hazards assumption. You can use the Kaplan-Meier estimator to get a non-parametric estimate of the discrete hazard (i.e. it does not make the proportional hazards assumption). Then you can log-transform those hazards to get a plot like the second. There are also statistical tests and other diagnostics (Schoenfeld residuals) to check this assumption.
The next question is do those violations matter? This somewhat depends on the question you are trying to answer and how substantial the violations are. In cases where the two log-hazard curves don't cross (like in D), it may be acceptable to smooth over the variations. However, it settings where the curves cross in the final figure, it is probably best to not ignore this problem.
The original write-up of this answer was here


How to reverse complement the DNA sequences for given inverse/reverse coordinates?

Since you tagged python, here is a solution that solves this in a single script using Biopython:

import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

Output on your example:

>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC


Interpreting large z scores from microarray data

This is a statistical question. What a higher z value means is that it is more extreme (if you assume a normal distribution), thus a more extreme value is less likely to happen by chance. Which is translated to a (lower) p-value.

What it is significant and not is arbitrary. Prior to the experiment you must decide which is your threshold, how much chance you allow. Apparently you decided to select 0.05 as your threshold, which means that one of 20 results might be due to chance. But your extreme z score is not more significant, it is less expected (under the null hypothesis, those p-values might change if you change your null hypothesis, ie you think your genes expression don't follow a normal distribution) but as significant as your 0.0499999 p-values.


You can simply replace the sequence and quality scores with a single * in each field. The CIGAR information will still be there, so it will be apparent where there were InDels and possibly mismatches (assuming your aligning uses = and X CIGAR operations), but since the actual sequence of these variants cannot be known this is quite anonymized.


How to compute the chance of failing to detect a gene given the detection limit of a protocol

A 90% loss can be rephrased as a 10% chance of detecting anything. So what we want to find is the probability of detecting 0 molecules, when we start with 7 and have 10% probability of success. Once can do that in R as follows:

> pbinom(0, 7, 0.1)
0.4782969

So ~50%, as they stated. I suspect that part of the confusion arises from the fact that the detection limit is due to lowly-expressed genes/transcript being heavily affected by this loss. So the probability of detecting a single molecule out of 4 original molecules (assuming 90% mRNA loss) is 34%, for 3 molecules it's 27%, for 2 it's 19% and for 1 it's 10%. I think the threshold of 5 is mentioned more because it's a nice round number than there's anything particularly different in detecting a gene with 4 vs. 5 molecules (34 vs. 41% probability).


finding RNA-protein physical interaction

If you were to treat the nucleotide like a ligand you could follow these steps on pymol, and find the measurements you desire.

Calculating distances to centroids is a little bit more complex, as you have to make a pseudo atom to measure against.


Short template DNA is a big problem with nanopore sequencing; it doesn't easily show up on a gel if it's well distributed across different lengths (as appears to be the case here from the read length plot), and sucks up both adapter and sequencing capacity.
The usual way to deal with it (for reads in the size ranges seen here) is to reduce the bead concentration when doing magnetic bead cleanup. Try something like 0.6X beads (e.g. 30μl of Ampure XP beads added to 50μl template DNA solution), which should still allow most of the long template sequences to bind, but leave shorter reads in solution.


What is the meaning of the "*" character in bwa fastmap's output?

It means the read is unmapped. No origin in the reference could be found that was sufficiently similar to the read to call it a proper alignment.

In the SAM specification, * is always used if the information for that field is not available. Therefore, if no mapping information is available (=unmapped) one sets *.


There is a pretty good description here. The position in the reads means literally the position in the reads, so position 2 means the quality of all bases on position 2 in the fastq file (all reads). In the description link, you can find how the boxplots and blue and red lines are calculated (and what they mean). They are all based on the total quality values of all bases on a certain position, all bases means all bases of all the reads on a particular position.


Select top 100 genes ranked by variance in read counts

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it m:

sel = order(apply(m, 1, var), decreasing=TRUE)[1:100]

sel then contains the indices into your matrix (or the original dataframe).

BTW, I hope you're not filtering by variance before performing differential expression. That would be incredibly incorrect. Doing this sort of filtering for performing PCA or something like that can be useful, though you should really be using normalized counts (or a transformation of them).


Smith-Waterman Pairwise Local Alignment Algorithm

There are many implementations of the Smith-Waterman algorithm. However, they are not all guaranteed to give the same output because although it's an optimal local alignment algorithm, it still depends on the chosen scoring scheme (i.e. substitution matrix and gap penalties).

With this in mind, you could try EMBOSS Water online:
https://www.ebi.ac.uk/Tools/psa/emboss_water/

Or in Biopython, the Align.PairwiseAligner() (which uses SW by default when aligner.mode is set to local):
https://biopython.org/DIST/docs/api/Bio.Align.PairwiseAligner-class.html


The general steps are:

Install deepTools
Run plotCoverage on your BAM files
Remove the left-most plot in photoshop/gimp/imagemagick


Converting a VCF into a FASTA given a reference with Python, R

You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the sequence, then output in FASTA format with whatever unique identifier for each sequence you want to use.

Edit: The FastaVariant class is new in pyfaidx, I hadn't seen it before. It makes this even easier, as you won't have to manually mess with the VCF to get the positions. Something like the following should work:

consensus = FastaVariant('yourgenome.fasta', 'variants.vcf.gz', het=True, hom=True)

out = open("variants.fasta", "w")

for chrom in consensus.keys:
    for var in consensus[chrom].variants_sites:
        record = consensus[chrom][var-1:var]
        print(record, file=out)

out.close()


How do I rename output folder created using WDL's glob function, with cromwell execution engine?

Ok, thanks to Ruben Vorderman on the openwdl github page, I found a solution that does not involve renaming the glob folder (which I believe is not possible?). The only way is to write out the names of the files and pass them explicitly to the Array[File] output. This is how I did it (simplified), in case it helps anyone else in the future:
task fastqc_task {
  File forwardReads
  File reverseReads
  String fwdName = sub(basename(forwardReads),".f.*q.*$","")
  String revName = sub(basename(reverseReads),".f.*q.*$","")
  command {
    mkdir -p ./qc/fastqc_zip
    fastqc ${forwardReads} ${reverseReads} --outdir ./qc/fastqc_zip
  }
  output {
    Array[File] fastqcArray = ["./qc/fastqc_zip/${fwdName}_fastqc.zip", "./qc/fastqc_zip/${revName}_fastqc.zip"]
  }
}


vcf file to gene annotations and or genes

Creating a FASTA file by applying variants found in a VCF file is called "creating a consensus FASTA file". You can use bcftools consensus. Unfortunately, the tools only works genome-wide (at the FASTA level), so you might need to use a different tool to extract a gene sequence from the FASTA file using the gene's co-ordinates.

If the VCF file is annotated with HGVS notations for the variants, there might be a way to use those notations to apply the edit to the reference gene/transcript sequence directly. I know the hgvs python module works well with the HGVS notation, but I don't see a way to "apply" the variant to the reference sequence.


Is there a "definitive" database for Mobile Genetic Elements?

Would advise doing Denovo identification method structure-based. Then use these public databases (  RepBase, Dfam ) to look for high similarity or homology exploration, this possibility if your species, the repeats not yet discovered, or deposited in the databases.

If the species well-investigated earlier and the mobile element recognised, I assume it will be a straightforward task using species-specific search via Repeatmasker tool with the -species flag to obtain the MobileElement annotation in the genome.

Let me know if you need any data from Repbase related to your species, i have an entree.


Is is possible to predict ncRNAs from sequence and homology alone?

Is there any conserved secondary structure in this RNA? One could align all known instances of this RNA and see if there is any significant covariation according to R-scape. That could serve as a computational signal that the RNA is real (although wet lab experiments would still be needed).
Another question is whether this RNA is already in Rfam? If so, then the job of making an alignment of known instances is easier.


Curl error when installing Bioconductor packages with Microsoft R Open 3.5.0

I can't comment yet but from what you said it seems that you need to install additional libraries like the one you mentioned: "curl".

from the R terminal you should type:

install.packages("curl")

then R will install the package. Not sure if there are some other dependencies. If yes R will return you some error message in which there are some package names you need to install.

EDIT:

of course, after installing all the required packages you should be able to install biocLite too.


I just copied your FASTA sequence into a file (on Mac OS) and uploaded it to Primer BLAST, the search did not give an error.

If the error is caused by the line endings, maybe you can fix your downloaded files with a tool like dos2unix.

Update:

On a whim, I tried adding a description to the FASTA header:

>test1 description
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAA
AATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGG
ACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTT
CGTCCGGGTGTGACCGAAAGGTAAGATGGAGAGCCTTGTCCCTGGTTTCAACGAGAAAACACACGTCCAACTCAGTTTG

This time, there was no error.

Does that work for you as well?


Retrieving Metabolic networks from WIT Database

It seems like there are some files related to the databases at the ftp site: http://ftp.mcs.anl.gov/pub/WIT2/. The related files were last updated on 2002.

I found it via searching at google: site:*.mcs.anl.gov/ WIT.

Which reports also a page about the database. There seems to be a project to develop WIT3, because there is a mailing list for developers, maybe you can ask them more details if you have problems with the database.


Which quality score encoding does PacBio use?

PacBio does use PHRED 33, but it turns out the question may be irrelevant for the newer PacBio Sequel Sequencer, because it reports all base qualities as PHRED 0 (ASCII !). The RS-II reports PHRED 33 quality scores.

The Sequel provides data in the form of subreads, which are the circular consensus sequences (CCS) from a single zero-mode waveguide (ZMW). Essentially, a single circularized DNA molecule enters a ZMW well, and is sequenced as many times as possible. (Sometimes multiple molecules get into a single well, but that's not the normal case.) The number of times this single molecule can be sequenced depends on several variables, including molecule length. The consensus sequence for this single molecule is then reported out as a subread.

As I started playing with the quality scores in my data, I realized the quality score for every base is !, which equates to a PHRED quality of 0. I contacted someone at PacBio who informed me that they learned the base quality values were not any more useful than the bases themselves for generating a high-quality circular consensus sequence. Thus, rather than spend the computational effort to calculate the scores, they skipped it, and reported ! for every base. They also report an overall read quality of 0.8, but this is simply a placeholder. They are working to make this information more readily available.


Ragout worked well for me. It is a reference-based scaffolder.

ragout.py --outdir output/ input/ragout.rcp > message.txt 2> error.txt --threads 8 --overwrite

The ragout.rcp file contains a 'recipe' with the reference sequence and the contigs in FASTA format.:

.references = ensembl
.target = assembly

ensembl.fasta = ref.fa
assembly.fasta = contigs.fa

In my case, it reduced the number of contigs 10-fold and increased the N50 10-fold. It added N's between the contigs.


Extracting p-values using corrr package in R

I would suggest looking at some of the documentation.
As far as I can tell, the intention of the package is to visualize and organize correlation coefficient estimates. It does not seem that the author of the package has any interest in p-value estimation.
If you are interested in computing p-values, I would suggest using a different tool, such as the builtin cor.test() function.


How hard is it to clean and QC gene expression microarray data?

I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There are many tutorials and review papers about normalization and QC of arrays. Try there and see if you can easily understand what these papers are telling you.

Read the limma user's guide, and try to follow the tutorials therein. If you understand that quickly you can try a real experiment of your own. If you don't understand what's described in there, you might want to look for another project for your thesis.


Parallelizing WGCNA k-means clustering and merging smaller clusters

The bad news is that, indeed, projectiveKMeans is not parallelized and I am not sure how much of it is (easily) parallelizable. The good news is that with 15k features (genes) and 96GB RAM you don't need to run preclustering at all. Just analyze the whole data set in one block.

If you for some reason insist on splitting the data into blocks, make sure the BLAS R uses on the cluster is a fast BLAS, ideally one that can parallelize matrix operations. A lot of time in merging of small clusters is spent calculating singular value decompositions, which becomes much faster using a good BLAS. HPC clusters will often have Intel KML available, or the admins could install OpenBLAS.


Indeed the web interface does not seem to link the full translation table, but it is available on their ftp server at releases -> release_ID -> precomputed_files -> genes -> fbgn_fbtr_fbpp_fb_2019_06.tsv.gz

For example, this is the latest (FB2019_06) translation table available.


panic: division of zero by zero or infinity by infinity error when running GODON positive selection analysis

Ah wait you mentioned a branch sites test, right that makes sense. You've got two or more identical sequences in the analysis. Ensure there are no identical sequences in your alignments and ... PRESTO! It will work :-)
Below are my original thoughts before I re-read your question. I've kept the notes because it helps explain the theory a bit. There is no need to supply the tree, it should be resolved.

The tree likelihood looks to be decreasing nicely, so it can see both the tree and the codons in the alignment.
At a guess the error relates to two or more indentical nucleotide sequences in your alignment. Thus zero non-synonymous sites (amino acid change) and zero synonymous (silent) sites. This error would only occur if you were calculating selection for each 'node' in the tree (internal branch).  Alternatively, possibly if the number of synonymous sites was zero, regardless of the non-synonymous mutations this error could appear, but this doesn't seem to fit the error output.
I don't know "Godon", but recommend PAML for this analysis. Admittedly it is complicated to deploy, but it is a very robust algorithm and hugely respected.


Devon will probably add a more precise answer, until then I just link his Tweet from last week:
https://twitter.com/dpryan79/status/1394714988720308226


Get a single number representing the contact probability between a pair of genomic loci using HiC data

The fundamental unit of a Hi-C contact is the nearest restriction enzyme cut site.
This is where each mate in a read pair from a Hi-C library comes from and is mapped to when converting sequencing data into a contact matrix.

Because the distance to the nearest cut site is variable at different genomic positions, and depth of coverage of any single locus is typically low, contacts are aggregated into bins and "normalized" to account for differences in cut site densities within those bins.
This normalization transforms the raw contact matrix into something resembling a stochastic transition matrix, where you can think of the elements of that matrix as "what is the probability that at any moment in time bin $i$ is in contact with bin $j$".

So the probability you're looking for is the $(i, j)$th element of the normalized contact matrix, where bins $i$ and $j$ contain your regions of interest.
If your regions straddle the bin boundaries, you may want to average those bins.

You get to pick the resolution of the contact matrix, but if you select a resolutiong that's too small, you run the risk of making your matrix too sparse and having no contact between those loci.
If you pick one that's too large, your matrix elements aren't really informative for your question at hand.

Given your 500-5000 bp range and median of 1 kbp, it may make sense to pick a resolution between 1 kbp and 5 kbp and use that.


How can I transform a mapped BAM file into an unmapped BAM file?

Looks like you need Picard's RevertSam operation (described in context here, manual here).

java -jar picard.jar RevertSam \
     I=input.bam \
     O=reverted.bam


How can I get unmatched reads for defective genomes analysis using bwa and samtools?

If you want to keep track of the unmapped reads, I recommend you to use Bowtie2 instead of bwa.
Start by building an index:
bowtie2-build genome.fa genome.fa

Then map your reads, using arguments to save unmapped reads in a specific file like this:
bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz

It will write the alignment as well as unmappedreads_1.fq.gz and unmappedreads_2.fq.gz that you may use in DI-tector.


Do any computational phylogenetic methods enable the specification of ancestral states?

For maximum parsimony, Fitch's method seems simple and intuitive, and can be done in the R phangorn package: Fitch, Walter M. (1971). "Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology". Systematic Zoology. 20 (4): 406. doi:10.2307/2412116. ISSN 0039-7989. JSTOR 2412116

There are various methodological criticisms of parsimony that you can find in the literature, but its appeal is you can easily understand why the algorithm came to the conclusions it did.

RAxML can do maximum likelihood reconstructions based on a tree and tip values, and I think Mr. Bayes can do a Bayesian version.


Download a release (Hammock_v_1.2.0.7z) from https://github.com/krejciadam/hammock/releases
Unzip
To run the first example:

cd /your_path/Hammock_v_1.2.0/examples/MUSI/
java -jar ../../dist/Hammock.jar full -i musi.fa

As mentioned in the comments, it's in the manual and any further issues need more information (errors etc).


Why sequence the human genome at 30x coverage?

The earliest mention of the 30x paradigm I could find is in the original Illumina whole-genome sequencing paper: Bentley, 2008. Specifically, in Figure 5, they show that most SNPs have been found, and that there are few uncovered/uncalled bases by the time you reach 30x:

These days, 30x is still a common standard, but large-scale germline sequencing projects are often pushing down closer to 25x and finding it adequate. Every group doing this seriously has done power calculations based on specifics of their machines and prep (things like error rates and read lengths matter!).

Cancer genomics is going in the other direction. When you have to contend with purity, ploidy, and subclonal populations, much more coverage than 30x is needed. Our group showed in this 2015 paper that even 300x whole-genome coverage of a tumor was likely missing real rare variants in a tumor.

On the whole, the sequence coverage you need really depends on what questions you're asking, and I'd recommend that anyone designing a sequencing experiment consult with both a sequencing expert and a statistician beforehand (and it's even better if those are the same person!)


How do population genetics people define a population?

A population is defined as a group of people (or organisms of the same species) living together in a geographical region. However, for human genetic studies, it is wise to use ethnicity as well (so group of people with the same ethnic background). See here how they define population in the 1000 genomes project.


RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'

This is resolved by installing the more recent version of the code from GitHub.

git clone https://github.com/lmcinnes/umap
cd umap
pip install --user -r requirements.txt
python setup.py install --user

As Peter said, this kind of question is best directed to the Issues page on their GitHub repository.


What methods should I use from PythonCyc API to query metabolites in BioCyc database?

UPD: As suggested by Llopis, maybe the correct identifier is VAL?

I think there is no metabolite with such name in the database. I tried their search, the correct name seems to be L-valine. I expect that they use similar naming for other amino acids. If you cannot find something, you can use their compound search through the website.


How to interpret UniProt allele patterns?

These are cytotoxic T-cell HLA alleles. HLA genotyping is very common and easy to do, so Genbank is the repository.

A*11:01 has very high frequency in Aborigine populations here . You can explore the population genetics and the past, present and current population genetics distributions per population of HLA at http://www.allelefrequencies.net/


Explanation for RNA-seq samples not clustering in PCA as expected

It's hard to say what is going on without looking at the plots. But it is a reasonable guess that samples are mislabelled and I think it is worth checking the experimental log carefully. It is also possible that the two samples are of low quality. People often use PCA to check and remove outliers from their RNA-seq libraries. So you may also check the RNA-seq library quality for those two samples.


Should biological replicates be the most similar pairs in an RNAseq experiment?

In summary you basically have noisy data and to simply answer your question it depends on the genetic diversity of the samples you are dealing with. Note, from my recollection this was confusing to start with because the "replicates" column was labeled as "treatment", but I could be wrong.

Revised answer, firstly there is zero stationarity - thus its a straight time-series.

Summary As a time series

The statistical noise between A and D is lower than the difference in phenotypic expression between each and all othe replicates in the analysis ABCD - FOR DAY 0 VERSUS DAY
4, so they cluster and are distinct.
The statistical noise between B and C is greater than the phenotypic variation between their expression profiles for DAY 0 to DAY 4. I suspect they are genetically very similar. Thats the conclusion, you can spice it up with different analysis, but I'd leave it because might find a lack of statistical significance in your core observations (A and D).
Within the time-points the data is very noisy and doesn't appear to show a distinct signal. Thats okay, this is essentially a time series.
Sample D is distinct in all scenarios (which is good)

Result To reiterate and expand a bit, the result you have is that sample A and D biological replicates are consisent within the analysis you are performing, i.e. day 0 verses day 4 (below). Thus they cluster together and this is proof of distinct phenotypes above the statistical noise (A with A and D with D). An intra-time point dendrogram would have a different result. The replicates samples B and C are not consistent and the interpretation is that the phenotypic expression profiles are so similar they are essentially identical within the boundaries of statistical noise.

Analysis The dendrogram is combines day 0 to day 4 analysis with day 0 to day 0 analysis and day 4 to day 4 analysis. HOWEVER, the clustering algorithm will correct for this and separate based on the largest distance (day 0 versus day 4). The signal that you get in your dendrogram is the largest signal in the data set which is day 0 vs day 4, and this is how UPGMA clustering works (I probably have previous posts on the statistical algorithm of UPGMA). So actually that is fine.

Advice With data this noisy I'd go no further in attempting to deploy critical probablity (a p value)  because this data could be the result of a complex distribution. This would require a statistician. Personally
 I'd upload the PCA as a separate post, but its just my opinion, because UPGMA versus eigenvalue based PCA are statistically quite different approaches.


Question about Co-expression analysis and finding targets for lncRNAs

You have several options to approach this with WGCNA (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.

Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).

I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.

I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.


You are printing the attribute before setting it, that's why it's complaining. So instead of this:

    print(i.attributes['Note']=["Gene description"])

Set it first and then print:

i.attributes['Note']=["Gene description"]
print(i)


Obtaining HGDP project data in fasta format

You can download HGDP data in FASTQ format here: https://www.internationalgenome.org/data-portal/data-collection/hgdp


Library for computing BWT-based alignments

First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.

It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).

Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).

There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:

SeqAn. See the FMIndex tutorial and the Pairwise Sequence Alignment tutorial for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a quick development of a read mapper, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in YARA.
SDSL-Lite. This is a general library for succinct data structures. See the tutorial slides for an idea of how it works. For instance, GramTools are based on SDSL.


Retrieve all genes for a specific gene ontology

You can use BioMart to filter by your GO term and get the genes as attributes. BioMart is ontology-aware so will pull out all genes associated with your term and with any of its child terms.

There is no need to look up the GO child terms as BioMart already deals with this. A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.

If you're not already familiar with BioMart, there's a help video.

Here's a BioMart query that gets you exactly what you need in human using the web interface

If you want to work with the R interface to biomaRt, here's some documentation and this is the identical query to above:

library(biomaRt)
ensembl=useEnsembl(biomart="ensembl", dataset='hsapiens_gene_ensembl')
getBM(attributes = c('ensembl_gene_id','external_gene_name', 'go_parent_term', 'name_1006'),
      filters = 'go',
      values = 'GO:0030098',
      mart = ensembl)


Preparing binary matrix data for Scikit classification algorithms

Since you are using R, you probably don't want to use scikit-learn, which is for Python. However, there is a similar R library mlr ("R package to make machine learning in R easy") that provides a unified interface to all popular machine learning methods. Check their tutorial on how to get started.


My current solution is (awk is not a tool I'm very familiar with, so this may not be good awk code, there is presumably a way to get it to print the header lines unmodified.):
samtools view -H temp.sam > temp2.sam
samtools view temp.sam |
awk '{OFS = "\t"} {for(i=1;i<=NF;i++)  if ($i ~ /^RX:Z:/) {$1=$1"_"$i; gsub("RX:Z:","",$1); print}}' >> temp2.sam

This seems to work when I test it, but I don't know if I'm missing anything.


What changed wwith Python in WDL change from draft 2 to version 1?

Complements of Phil Shapiro;
Use ~{} instead of ${}.
So this example from the current docs is wrong:
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("${in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}

It should be :
task heredoc {
  input {
    File in
  }

  command<<<
  python <<CODE
    with open("~{in}") as fp:
      for line in fp:
        if not line.startswith('#'):
          print(line.strip())
  CODE
  >>>
}


To add to the answer from Bioathlete, "reads" in the traditional genomic sense are the string of nucleotide data inferred from the sequencer from a given fragment of genetic material (could be Illumina, Nanopore, PacBio, 10X, etc.).

So read count first applies for all of the "reads" defined above generated by the sequencer. From a NovaSeq, the total number of bases from a single flowcell is up to 3 trillion bases at 150bp each read for the highest throughput flowcell (so a read count of 3,000,000,000,000/150 = 20B total read count).

The read count for a single Nanopore flowcell is roughly (30 billion bases/avg read length), but the read lengths can be SUPER long (100kb compared to 150bp is substantial). Long read technologies are generally used for different analyses than short reads.

So all of the above reads would be defined as "raw reads", meaning they are reads not yet analysed in any way to be used for data analysis. So using illumina as an example, for a single Whole Genome sequencing sample we might have 800,000,000 raw reads. We then use an aligner like bwa mem to map the reads to a single reference genome (or in your case you would use a different aligner to map to multiple genomes). From here we still have the concept of raw reads, but we also now have the total number of mapped reads. A typical good illumina library for a single reference genome will have >90% mapping rate for a human sample.

We can then say we have 800,000,000 raw reads, and 720,000,000 mapped reads with a mapping to a reference genome of 90%. So we can compare "raw read count" to "mapped read count" to get an idea of the quality of our library.

In the case of metagenomics, you would consider your "raw read count" to still be the total number of reads you got off the sequencer, but the "mapped read count" would be the sum of the reads mapped to each individual reference genome to get a total quantity of "mapped reads" which could give you an indication of library quality. I'm not sure what a good metagenomics library "mapped read count" is, but I'm sure you can find that information.

Lastly, we want to get an idea of the depth of coverage for each mapped read to the genome. This will tell us how well our mapped reads actually covered the reference genome or genomes we are interested in.

I hope this was helpful!


Can I submit a R package to Bioconductor or CRAN if I have already published it a journal?

You can submit a package to a repo after the publication, but I would say that you really should do it before.

Bioconductor accept only packages that are not published on CRAN, however, academic publications are fine.

Conversely, Bioinformatics does not enforce any particular platform for sharing the code, but you have to make it available somehow (e.g. GitHub repo, university webpage, bioconductor). From the instructions for authors

  If the manuscript describes new software tools or the implementation
  of novel algorithms the software must be freely available to
  non-commercial users at the time of submission, and appropriate test
  data should be made available.

However, you as a developer should make all the effort to make your software available and easily accessible. Also, sooner you make it available, sooner people start to cite you, which is a great advantage for you. Finally, you really need users to make sure that your package really works (i.e. debugging all the corner cases).

I could not emphasise more, the correct order is 1. make it available 2. write a publication


It looks as though the script is looking for a URL (fastaURL) so you could try specifying the file as file://directory/test.fasta

By default though, javascript cannot read local files from disk for security reasons. Are you sure it is even possible to specify a local file in your script in this way?  If yes, then the problem is likely that you need to specify the path correctly.


Why use Needleman-Wunsch if there is no way to evaluate the statistical significance

It's best to flip the question on its head and ask yourself instead, "Do we care about the statistical significance of an individual alignment?" The answer is "yes" if you're in a context like homology modeling, where you're aligning two protein sequences. The answer, however, is a very definite "no" in the context of NGS, where you just want the most likely alignment and a vague assessment of its likelihood of being wrong (for filtering purposes).


Create a NHX phylogeny format with r or python

The treedata book has a section on tree annotations that should work for you.


There's no equivalent to the wiggle header in bigWig (or bigBed) files, which is why UCSC uses the file name. This is actually the reason for the track line stuff that you linked to, since you can then specify a name and just point to where the bigWig (or other format) file is on the internet.

BTW, you can certainly convert your bigWig to wiggle, add the track line, and convert it back with UCSC tools (e.g., bigWigToWig and wigToBigWig). It shouldn't change the track name when you display it in UCSC (though if it does, then let me know, since I'll modify pyBigWig and libBigWig to support that then).


RPGC normalisation creates artefacts at centromere

UPDATE

Here is some other protein ChIPs that I ran through the same normalisation:

As you can see, they don't seem to be enriched at the centromere. So unless anyone brings other contributions, I will accept this as answer.

RPGC normalisation doesn't appear to create artefacts at centromere. The peaks observed here are as real as can be.


are GSEA and other geneset enrichment analysis supposed to yield extremely different results between them?

Each of these methods do something different, so it is reasonable to expect different results. The bottom line is that there isn't a single question you can ask when you do an "enrichment test". For example: you can test if a set of genes are more expressed than this other list of genes or they are above the mean levels of any other gene.
Even more, you use the gene ontology, which has an underlying structure (directed acyclic graph) that no method you used takes into account but could modify the results if taken into account. However, the KEGG pathways do not have this structure. So when you want to apply a method mixing both data sources you need to decide if you want to use it or not.
This is also complicated because for most of these methods you can only compare the p-value which depends on the underlying assumptions of the methods (uniform, normal, binomial distribution?)
The latest comprehensive study comparing enrichment methods including some of the methods you used and more is this one (that I know of): Toward a gold standard for benchmarking gene set enrichment analysis.
Last, what would a correlation tell you here? It is know that there is little agreement between methods. You couldn't say that one method is equal to the other unless you tested with several datasets and you had always a correspondence on the p-value.


VCF-merge fails due to tabix not producing .tbi files

You should go with the CSI index as it seems like the BAI indexes cannot handle chromosomes longer than 2^29.

Your problem is that you're trying to merge the index files using vcf-merge. You should merge the VCF files, not the .csi index files.


How to transfer gff annotations in genome with extensive duplications?

There is one very simplistic way I use which might work for what you are doing, it is similar to what terdon proposed.

Take a de-novo microbial genome annotation tool (I have my own, but you could use/modify prokka). Tools like these often first predict gene boundaries (with other tools like prodigal or glimmer) and then try to assign a function to found genes. This function assignment is often done with BLAST and other tools ... and that is where you can go in and modify to do what you need.

I use a "knowledge" protein database of genes I want to have very strictly annotated as a first line of annotation (e.g. in your case: the annotated genomes). For that I loop through very strict identity/similarity parameters which get gradually relaxed.

E.g.:
Loop 0: only transfer annotations at 100% DNA identity, same length.
Loop 1: only transfer annotations at 100% similarity, same length.
Loop 2: only transfer annotations at 99% similarity, length +/- 1%.
...
Loop n: only transfer annotations at 100-(n-1)% similarity, length +/- (n-1)%.

In each loop, obviously only annotate what has not been annotated in previous loops.

After that, use "normal" annotation pipeline of the tool to annotate the rest.


What are the largest fully publicly available molecular QTL datasets?

The largest eQTL datasets I can think of is probably eQTLGen. This is a meta-analysis of many studies using >30,000 samples from whole blood. This might be less ideal if you want cell-type-specific QTL information.

For other molecular QTLs, I do not think there are large-scale datasets yet.


Database of position weight matrices for protein motifs?

PFam has HMM models for protein domains, this is almost certainly what you're after. You can download all PFam-A models from ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz.

Then, you're looking for the "Pkinase", PF00069 family (I found this simply by searching the PFam website for "kinase"), so extract it from the file:

zcat Pfam-A.hmm.gz | awk -vRS='//' -vORS='//' '/PF00069/' > PF00069.pfam

Then, format it for HMMER format (you will need to install HMMER)

hmmpress PF00069.hmm

And now you can use it to search your sequences (fasta or flat file, probably others; see the documentation for details):

hmmscan PF00069.hmm prots.pep

I tried this on a fasta file with the sequences of CBK1_YEAST, CSK21_YEAST (kinases) and FLNA_HUMAN (non-kinase control) and got:

# hmmscan :: search sequence(s) against a profile database
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:             prots.pep
# target HMM database:             PF00069.hmm
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       sp|P53894|CBK1_YEAST  [L=756]
Description: Serine/threonine-protein kinase CBK1 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CBK1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    2.7e-64  203.2   9.0    3.5e-46  143.8   0.2    3.1  3  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ?   -2.9   3.2      0.19      0.19     198     226 ..     222     243 ..     191     277 .. 0.57
   2 !  143.8   0.2   3.5e-46   3.5e-46       1     153 [.     352     505 ..     352     522 .. 0.95
   3 !   63.6   0.0   9.8e-22   9.8e-22     154     264 .]     568     672 ..     551     672 .. 0.83

  Alignments for each domain:
  == domain 1  score: -2.9 bits;  conditional E-value: 0.19
                           TTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
               Pkinase 198 sgekgkekvekeldqlekilkilgetkek 226
                           + +++++++++      + ++i ++++++
  sp|P53894|CBK1_YEAST 222 Q-QQQQQQQQQ------QHMQIQQQQQQQ 243
                           1.111111111......122222221222 PP

  == domain 2  score: 143.8 bits;  conditional E-value: 3.5e-46
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHH.HHHHHHHHHHHHHSSTTB--EEEEEEESSEEEEEEE--TTEBHHHH CS
               Pkinase   1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkke.kkvlrEikilkklkhpnivrlyevfeekdelylvleyveggeladl 88
                           +++++ +G+G+fG+V  + +k+tgki+A+K++ k+++ kk++  +v +E  +l+ +++p +v ly++f++ ++lyl++e+++gg+l+++
  sp|P53894|CBK1_YEAST 352 FHTVKVIGKGAFGEVRLVQKKDTGKIYAMKTLLKSEMYKKDQlAHVKAERDVLAGSDSPWVVSLYYSFQDAQYLYLIMEFLPGGDLMTM 440
                           67899********************************999999********************************************** PP

                           HHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEETTTEEEE-SGTTSEECCSSCS CS
               Pkinase  89 lsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLidekgelKitDFGlakelesssk 153
                           l r + ++e+ +++++++ + ++e +H+ g+iHrD+Kp+NiLid +g++K++DFGl++ ++++++
  sp|P53894|CBK1_YEAST 441 LIRWQLFTEDVTRFYMAECILAIETIHKLGFIHRDIKPDNILIDIRGHIKLSDFGLSTGFHKTHD 505
                           **999*****************************************************9999876 PP

  == domain 3  score: 63.6 bits;  conditional E-value: 9.8e-22
                           BCTCCSCGGGS-HHHHTTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH......TTS-HHHHHHHH CS
               Pkinase 154 ltsfvgtreYlAPEvlkeneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkeklpeaselseeakdllk 242
                           + s vgt++Y+APE+  +++y++++D+WslG i+ye l g ppf +e+ +e+++k+++          e+++++p + ++s ea dl++
  sp|P53894|CBK1_YEAST 568 AYSTVGTPDYIAPEIFLYQGYGQECDWWSLGAIMYECLIGWPPFCSETPQETYRKIMNF---------EQTLQFPDDIHISYEAEDLIR 647
                           6789*****************************************66665533333333.........557799*************** PP

                           HHT-SSGCCSTT....HHHHHTSGGG CS
               Pkinase 243 kllkkdpkkRlt....aeellqhpyl 264
                           +ll+    +Rl+    a+e+++hp++
  sp|P53894|CBK1_YEAST 648 RLLTHA-DQRLGrhggADEIKSHPFF 672
                           ***875.678877779********98 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (756 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.00
# Mc/sec: inf
//
Query:       sp|P15790|CSK21_YEAST  [L=372]
Description: Casein kinase II subunit alpha OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CKA1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    7.5e-62  195.1   7.2    3.2e-55  173.4   0.7    2.1  2  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 !   25.5   1.1   4.3e-10   4.3e-10       1      60 [.      40      94 ..      40     115 .. 0.87
   2 !  173.4   0.7   3.2e-55   3.2e-55      56     264 .]     128     363 ..     101     363 .. 0.92

  Alignments for each domain:
  == domain 1  score: 25.5 bits;  conditional E-value: 4.3e-10
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHHHHHHHHHHHHHHHSSTTB CS
                Pkinase  1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkkekkvlrEikilkklkhpni 60
                           ye ++k+G+G++++V+++v+ +++ ++++K++k  k+kk k     rEikil  l+ +++
  sp|P15790|CSK21_YEAST 40 YEIENKVGRGKYSEVFQGVKLDSKVKIVIKMLKPVKKKKIK-----REIKILTDLSNEKV 94
                           78999****************************99888855.....7******9998877 PP

  == domain 2  score: 173.4 bits;  conditional E-value: 3.2e-55
                            S.STTB--EEEEEEE..SSEEEEEEE--TTEBHHHHHHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEE.TTTEEEE CS
                Pkinase  56 k.hpnivrlyevfee..kdelylvleyveggeladllsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLid.ekgelKi 139
                            + h ni++l++++++  +++  lv+eyv++ ++  l+    +l++ e+++++ ++l++l+y+Hs+gi+HrD+Kp+N++id ++++l++
  sp|P15790|CSK21_YEAST 128 NgHANIIHLFDIIKDpiSKTPALVFEYVDNVDFRILYP---KLTDLEIRFYMFELLKALDYCHSMGIMHRDVKPHNVMIDhKNKKLRL 212
                            35*************8878888*********9988887...5**************************************66677*** PP

                            -SGTTSEECCSSCSBCTCCSCGGGS-HHHH.TTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
                Pkinase 140 tDFGlakelessskltsfvgtreYlAPEvl.keneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkek 226
                            +D+Gla+ ++ + +++ +v++r +  PE+l ++++y+ + D+Ws G +l++++++++pf   +g++    ++dql+ki+k+lg+++++
  sp|P15790|CSK21_YEAST 213 IDWGLAEFYHVNMEYNVRVASRFFKGPELLvDYRMYDYSLDLWSFGTMLASMIFKREPFF--HGTS----NTDQLVKIVKVLGTSDFE 294
                            ******************************999***************************..3443....789**********77766 PP

                            ....TTS...............................-HHHHHHHHHHT-SSGCCSTTHHHHHTSGGG CS
                Pkinase 227 lpeasel...............................seeakdllkkllkkdpkkRltaeellqhpyl 264
                                +                                 ++e++dl+++ll++d+++Rlta+e++ hp++
  sp|P15790|CSK21_YEAST 295 KYLLKYEitlprefydmdqyirkpwhrfindgnkhlsgNDEIIDLIDNLLRYDHQERLTAKEAMGHPWF 363
                            5554444588999999****************************************************9 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (372 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.02
# Mc/sec: 4.91
//
Query:       sp|P21333|FLNA_HUMAN  [L=2647]
Description: Filamin-A OS=Homo sapiens GN=FLNA PE=1 SV=4
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------

   [No hits detected that satisfy reporting thresholds]

Domain annotation for each model (and alignments):

   [No targets detected that satisfy reporting thresholds]

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (2647 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         0  (0); expected 0.0 (0.02)
Passed bias filter:                        0  (0); expected 0.0 (0.02)
Passed Vit filter:                         0  (0); expected 0.0 (0.001)
Passed Fwd filter:                         0  (0); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               0  [number of targets reported over threshold]
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
# Mc/sec: inf
//

HMMER is a great tool and can also analyse multiple sequence alignments to build the hmms if that's what you need. I strongly urge you to reads its documentation which is both informative and a surpisingly fun read (it includes such phrases as "Their format is “proprietary”, which is an open source term of art that means both “I haven’t found time
to document them yet” and “I still might decide to change them arbitrarily without telling you”").


Why do I get so many insertions from Minimap2 on my Nanopore WGS?

Insertions and deletions are the dominant error mode of long read sequencing, including nanopore sequencing. What you see is not unexpected. Things may have improved by now if you would download the raw fast5 data and repeat the basecalling. There is no need to gunzip the fastq.gz prior to alignment. Your commands for alignment look alright, except that (if you are using a recent version of samtools) you don't really need samtools view and can pipe directly to samtools sort, using -o to specify the output format:

minimap2 -ax map-ont -t 20 <index> <fastq> | samtools sort -@5 -o alignment.bam

In a recent version of igv you can right-click the alignment and chose a "Quick consensus mode" (or something like that), which will make your alignments at least a bit prettier. Note how the consensus sequence (top of the graph) is quite clean - while individual reads may be noisy this is mostly randomly distributed and as such ironed out in consensus.


WGCNA co-expression network analysis with less than 20 samples

It's hard to say without knowing how different the subtypes are. If you do a common WGCNA, you may find modules related to the differences between A and B as well as to drivers of expression variation in A and B (common and separate). Since A has fewer samples, the variation will have to be stronger to be reflected by WGCNA modules.

You can run the WGCNA, find modules that vary the most in A samples (these may or may not also vary in B samples), then see whether LncRNAs belong to these modules and have a high kME (module membership).


How to find/build the evolutionary history of a protein from its sequence?

The general procedure is:

Find the sequence for the same (or as similar as you can find (see mentions of blast in the comments above)) gene in other species. You can use uniprot or any other large sequence database for this.
You now have a large set of sequences related to yours, which you'll need to compare. Search the methods section of your favorite paper to see how they reconstructed likely evolutionary history from this, there are a number of packages out there.


What is the difference between ENCODE Tier 1, 2 and 3 cell types?

The tiers denote the original priority of sequencing and processing the samples. This is understandable given the number of marks and experiment types that ENCODE tried to sequence. Further details are available on ENCODE's website.

Edit: Just to expand a bit, tier one cells were supposed to get processed first. Then tier 2 came second. Of course tier 2 had a LOT more cells, so it got split into 2 and 2.5, since there were already apparently tier 3 cells.


How to download sequence alignments against UniProtKB sequence database for a particular Pfam Id using R?

The following code served my purpose.
library(readr)
pfam_url <-
  'https://pfam.xfam.org/family/PF02171/alignment/uniprot/format?format=pfam&alnType=uniprot&order=t&case=l&gaps=default&download=1'
seqs <-
  read_table(
    pfam_url,
    col_names = FALSE,
    col_types = cols(X1 = col_character(), X2 = col_character())
  )


The scran::quickCluster method has a BPPARAM (i.e., BiocParallelParam), so if one provides something like

clusters <- quickCluster(sce, BPPARAM=BiocParallel::MulticoreParam())

it will use all available cores (a specific number of cores can also be given). See the intro doc to BiocParallel for info on other options.


Simple answer: Yes, that can very well happen depending on the average fragment size, the read length and the length of the UTR.


Yes, you can safely concatenate the technical replicates. Odds are good that these are even the same libraries just sequenced twice, so even labeling them as replicates is a bit of a stretch. As an aside, it would be surprising if you actually had a batch effect in a situation like this. You will commonly see sequencing facilities just sequence a given sample a second time if it didn't get the number of reads requested by the client to begin with. They then just have the client concatenate the files.


Your code looks fine, but rbind to progressively build and resize a data frame can be gluttonous in terms of memory and time, especially given the A-row x B-row times you are doing this operation, which is fine for very small datasets, but not often realistic for real-world data.

Instead, perhaps cat results to a text file, where results are written in tab-delimited columns. You can then read in that text file later on.

In addition to serializing or saving your results to text, which make them easier to work with in non-R contexts, you won't have to deal with R's memory allocation problems when running rbind repeatedly. You just write out the answer to a file and read it in later on, in one pass.

Another approach is to preallocate memory to a presized data frame, and then write to that allocated space.

You know the size of your inputs ahead of time, which makes this approach tenable. In this case, it looks like you are storing two floats and a string (rho, p.val, and Gene) and that you are going to store A-row x B-row number of them.

Here's one way you might go about this:

# create a presized data frame
Arows <- nrow(A.data)
Brows <- nrow(B.data)
correlation.result <- data.frame(t(matrix(c(0.0, 0.0, "foo"), nrow=3, ncol=Arows * Brows)), stringsAsFactors=FALSE)
colnames(correlation.result) <- c("rho", "p.val", "Gene")

# assign values to that presized data frame
Acols <- ncol(A.data)
Bcols <- ncol(B.data)
for (i in 1:Arows) {
  # get A data
  Ai <- as.numeric(A.data[i, 2:Acols])
  Aname <- A.data[i, 1]
  for (j in 1:Brows) {
    # calculate row offset in presized data frame
    ij <- (i - 1) * Brows + j
    # get B data
    Bj <- as.numeric(B.data[j, 2:Bcols])
    Bname <- B.data[j, 1]
    # run test
    correln <- cor.test(Ai, Bj, method = "spearman", exact = FALSE)
    rho <- correln$estimate
    p.val <- correln$p.value
    Gene <- paste(Aname, Bname, sep="_", collapse=" ")
    # assign value to row of result data frame
    correlation.result[ij,] <- c(rho, p.val, Gene)
  }
}

This avoids R needing to resize and recopy the data frame contents as they grow, which I suspect is the major reason it is taking so long to run.

Chapter or "Circle" 2 of Patrick Burns' book The R Inferno discusses the problems with growing objects with rbind and ways to handle it — you might find this a useful read. This book is a PDF download via your search engine of choice.


Does the DNA or RNA of some ethnic groups map better than others' to the human reference genome sequence?

For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.

In studies where they compare different ethnic groups (e.g., 1000 G) they rather not discuss such matters, but act like the analysis is 100% perfect.

There is a nice comment in Genome Biology about this issue.


how to change the PC use in the dimplot and feature plot

You can choose which PC dims to plot by specifying the dims argument within the functions. To use PC-4 and PC-5,

DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "pca")

and

FeaturePlot(object = pbmc_small, features = "Xist", dims = c(4, 5), reduction = "pca")

Seurat documentation is available here.


Subset smaller BAM to contain several thousand rows from multiple chromosomes

If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:

samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

That will select 40% (the .4 part) of the reads (123 is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the .4 part to a sufficiently small number.

In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.


Any fast options to query large VCF bed intervals?

this:

bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.


Difference between CPM and TPM and which one for downstream analysis?

You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-normalized values of the other genes).
If one has to choose between those two choices one typically chooses TPM for most things, since generally the length normalization is handy. Realistically, you probably want log(TPM) since otherwise noise in your most highly expressed genes dominates over small expression signals.


Finding peaks and estimating cell population sizes in multi-dimensional flow cytometry data

There are quite a few algorithms developed for the automatic classification of multidimensional flow cytometry data, you can see a (not so recent?) review here: http://www.nature.com/nmeth/journal/v10/n3/full/nmeth.2365.html. In your case, you are interested in the unsupervised methods, since you do not have data from single populations (which would provide you "labels" associated to each sample), as I understood.

In other words, if you have your data set as a matrix of 5 million samples (cells) with 14 dimensions (FC values) each, then you are looking basically for algorithms that will do a clustering of the samples into groups. By the way your description of each sample belonging to a multivariate Normal distribution leads to a Gaussian Mixture Model, which can be estimated from your data using e.g. the scikit-learn library for python: http://scikit-learn.org/stable/modules/mixture.html


I strongly suggest that you not try to come up with your own package for this when things like CSAW already exist in bioconductor and provide a number of useful normalization options.

For visualization purposes, I think it's best to simply take the log2 ratio (e.g., using bamCompare from deepTools) or alternatively normalize by total coverage (possibly excluding blacklist regions, if applicable for your organism).


creating population AF from set of individuals in 1000 genomes population?

bcftools has a plugin called fill-tags which should do the job for you easily, without having to mess around with lots of complicated bash scripting.

You could do something along the lines of this, presuming you have a .vcf file with all the IBS individuals in it:

bcftools +fill-tags in.vcf.gz -Oz -o out.vcf.gz

This will give each position in the .vcf file new INFO tags, one of which will be allele-frequency.

You can then pull off the frequency for a particular variant (CHR:POS) by:

bcftools query -t CHR:POS -f'%CHROM\t%POS\t%AF\n' out.vcf.gz