I normalized a high throughput dataset for a school project using DESeq library using the script bellow. The code is based on a lesson I had. My goal was determine the over expressed genes, but the normalization step should be the same.

#####################################################################
### DATA PRE-PROCESSING: Normalising gene expression distributions
#####################################################################
# Genes that are not expressed at a biologically meaningful level
# in any condition should be discarded to reduce the subset of genes
# to those that are of interest, and to reduce the number of tests carried out
# downstream when looking at differential expression.
# Using a nominal CPM value of 1 (which is equivalent to a log-CPM value of 0
# a gene is deemed to be expressed in a given sample if its transformed count
# is above this threshold, and unexpressed otherwise.
# Genes must be expressed in at least one group (or roughly any three samples
# across the entire experiment) to be kept for downstream analysis.

library("DESeq")

keep.exprs <- rowSums(cpm>1)>=3 # is filtering is ok
DGEfiltered <- DGEobject[keep.exprs,, keep.lib.sizes = FALSE]

## Normalisation generally required to ensure that expression distributions
## of each sample are similar across the entire experiment.
## Any plot showing the per sample expression distributions, such as a density or boxplot,
## is useful in determining whether any samples are dissimilar to others.
## Normalisation by the method of trimmed mean of M-values (TMM)
## is performed using the calcNormFactors function in edgeR.
# cpmF = filtered subset # log.cpmF = log cpm filtered subset
# Copy unnormalized, filtered data to show in boxplot
DGEf.notnorm <- DGEfiltered # normaliseerde
DGEf.norm <- calcNormFactors(DGEfiltered, method = "TMM")
DGEf.norm$samples$norm.factors

# --> scaling factors relatively close to 1 --> mild TMM normalisation

## Boxplots of log-CPM values showing expression distributions
## for unnormalised data (A) and normalised data (B) for each sample
par(mfrow = c(1,1))
par(mar = c(7,4,2,1))
# Boxplot unnormalised data
DGEf.notnorm$samples$norm.factors
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
log.cpmF.notnorm <- cpm(DGEf.notnorm, log = TRUE)
boxplot(log.cpmF.notnorm, las = 2, col = sample.col, cex = 0.9,
        main = "A. Unnormalised data", ylab = "Log-cpm")
# Boxplot normalised data
DGEf.norm$samples$norm.factors

log.cpmF.norm <- cpm(DGEf.norm, log = TRUE)
boxplot(log.cpmF.norm, las = 2, col = sample.col, cex = 0.9,
        main = "B. Normalised data", ylab = "Log-cpm")
# --> only when scaling factors are not close to 1 you will see a big difference


After completing the evaluation, how to predict PTM lysine or any other modification with unlabeled data for making a web server?

Interesting application. I would suggest using Python's Flask system, which is described in a few pages in O'Reilly's Programming PyTorch, chapter 8. PyTorch in Production. Django is not really recommended (don't say that on Stackoverflow however). O'Reilly do an entire book on Flask, but thats to make the site pretty.

Obviously, the explanation here is Deep Learning rather than ML but the basic gist will hold. You should also think about building a docker for distribution which is easy (same chapter). We're in lockdown so I don't have this book with me.

The difficult part will be web architecture because you shouldn't have the CPU/GPU on the same server as the quiry interface, otherwise you are prone to a denial of service attack or just sytem overload. To do this you need to consider  Representational state transfer (RESTful) and that is not easy at all (at least for me). Flask facilitates RESTful, but thats like systems admin stuff. I guess before getting into RESTful I'd see how popular your website is, but if its the same server as your institute website you'd need to be really careful because an external user could bring your institutional website down (you wouldn't be popular).


I don't program in R, but replace the last line of you code with the following and I suspect you will find it works,

for(i in seq(1,10000,1)){ recs <- entrez_summary(db="pubmed", web_history=pubmed_search$web_history, retmax=100, retstart=i) cat(seq_start+99, "sequences downloaded\r") }

Rationale I would guess seq is similar to the linspace function in Python, and put simply 10000 (total sample size) / 100 (the step size) = 100, 100 being the number of downloads you observed. In the code above its 10000/1. I would also suspect "step size = 1" is default, if this was correct then seq(1,10000) would give you your required output.

As a random aside I do like the presence of curly brackets defining the loop, it reminds me of Perl.


Gene name

Using an Ensembl ID, Uniprot ID or NCBI ID it is much easier to look up your gene. So I would strongly reccomend using one of those as a key. Personally, I'd use Uniprot, but Ensembl is probably better if your genes are human.
One reason for this is also that different resources will use different names, especially for humans.

There are a few considerations, so I will first say how this could be done properly...

Custom Uniprot matcher with weights

Some names are used by different protein, so which one is which would depend on context, co-occurrence with another gene.

Your dataset might be extremely noisy. Therefore, you likely have assigned a confidence score to your matches, yes? It would be nice to give a score to the name based on whether it is the first name of the protein or the fourth synonym in the list.

Therefore to do a conversion properly you really ought to make your own conversion system based upon context and Uniprot data. Uniprot contains the most data on synonyms and notes and binding partners etc. The dataset/API is in XML and Biopython does not have a parser, but I have some notes on the matter and can share some code if needed as I have done something similar (e.g. JSON for human gene/protein names to Uniprot IDs).

Lazy match

The NCBI entrez API is rather flexible and good at presenting to you matches for a strange names. So I would suggest using that.


OK I found the bug. CIBERSORT seems to require tsv files instead of csv files. Plus it doesn't seem to like empty columns


calculating nucleotide frequency per column


Error in checkFullRank(modelMatrix) deseq2

See the second example from the vignette, which is exactly your situation. The software cannot differentiate between the batch effects vs. condition effects because the groups are exactly the same.


Why is ribosomal RNA difficult to remove even with Poly(A) selection?

We've found ribosomal RNA to be less of a problem with sequencing that depends on polyA, which suggests the issue might be in the library preparation, rather than the selection.

Many polyA RNA library preparation methods involve amplification, rather than selection, which means that existing transcripts that are present in very high abundance (such as rRNA) will still show up in the mixture.

But polyA-tailed RNA isn't the only RNA present in a sample. For a comprehensive look at RNA, the whole-RNA sample preparation methods can't use any form of selection to identify target molecules. In this case, it's typically recommended to use a form of ribosomal RNA removal, where a probe to ribosomal RNA is attached to magnetic beads (or a column, or similar), and used to fish out out ribosomal RNA, leaving the remainder behind. Unfortunately those probes don't work perfectly, either due to the ribosomal RNA being damaged at the probe site, or the probe not matching the ribosomal RNA molecule, or the RNA just not getting into the right physical state to bind with the probe. All these things mean incomplete removal of ribosomal RNA, ending up with a situation where the highly-abundant ribosomal RNA still pokes through in the reads.


can we show the orig.ident in separate windows by using Seurat

The DimPlot() function of the new version of Seurat, Seurat v3 has a split_by argument, which splits the plot based on the levels of the variable provided. You basically need to provide a column from the @metadata slot.

If you are using Seurat v2, you can have the DimPlot() function return the ggplot object created and then use facet_wratp() or facet_grid() to generate facets based on your variable of interest.


Rapid metagenomics classifiers on long read data

You should also keep in mind that single read accuracy of ONT still is a bit lacking, a 2D or 1D^2 accuracy of 95% still means that on average there's one error every 20 bp, and due to the nature of the data some stretches may be more junky than others. Maybe centrifuge doesn't like that.

Getting a feel of what the data can tell you (and what not) is pretty important.

You may want to try other classification tools to get an idea how well centrifuge works, Kraken comes to mind. But also: what about BLASTing some of these reads at NCBI and simply look at the results to get a first impression?


It's not an unreasonable approach. rMATs is rather picky about its input, but it seems you've noticed that.
rMATs can't handle trimmed reads. It also can't handle soft-clipped reads. You might as well not trim. You'll get a lower alignment rate, but your only other choice would be to either exclude the trimmed reads or trim everything down to the same length.
You have a mixture of RF and FR alignments and the plotting script can't handle that.
Assuming you mean the read length distribution, it's actually shifted to the left, not the right. You had a single length before and a long left tail after trimming.


How to make a robust amino acid phylogeny with adequate robustness support

To finish you need bootstrap support using a model-based tree building algorhithm, I use RAxML, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,

raxmlHPC -b 459 -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

Make a consensus tree of the bootstraps,

raxmlHPC -f b -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)

You need a fairly robust desktop, dedicated for around 24 hours to complete one calculation. You have enough results to draw a conclusion.

Viewing the tree, FigTree is easy.

Root

Check the work of Cavalier-Smith, he heavily focused on basal protists to find a species which evolved, using external evidence, before any other member of your group. One of the ingroups may have sufficient evidence to say it evolved before the rest and that is your root, which you select as an outgroup.


Gene information from chromosome location

How can I quickly find out if it lies in an intron or exon of a gene?

I have an answer on Biostars that might help: https://www.biostars.org/p/124515/#124522

  Or it is in the upstream or downstream regions of a gene?

One option is to build a table of annotated bins and then map your regions-of-interest to those annotations.

These annotations would consist of whether the bin is inside a gene, or, after picking the nearest gene by strand, locate the bin to an upstream region of the gene TSS or to an downstream region of the gene TES.

To start, one could download stranded gene annotations:

$ GENCODE_URL=ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.basic.annotation.gff3.gz
$ GENCODE_GENES=${PWD}/../data/gencode.v28.genes.bed
$ wget -qO- "${GENCODE_URL}" | gunzip -c | awk '($3=="gene")' | gff2bed - | awk -v FS="\t" -v OFS="\t" '{ n=split($10,a,";"); p=a[4]; n=split(p,a,"="); $4=a[2]; print $0; }' | cut -f1-6 > ${GENCODE_GENES}

Chop up nuclear chromosomes of the genome of interest into 1kb bins (or whatever gradation of bin size you prefer):

$ ASSEMBLY=hg38
$ GENOME=${PWD}/../data/${ASSEMBLY}.bed
$ CHOP_SIZE=1000
$ CHOPPED_GENOME=${PWD}/../data/chopped.${ASSEMBLY}.bed
$ fetchChromSizes ${ASSEMBLY} | awk -v FS="\t" -v OFS="\t" '{ print $1, "0", $2 }' | grep -ve "\.*_\.*" | grep -v "chrM" | sort-bed - > ${GENOME}
$ bedops --chop ${CHOP_SIZE} ${GENOME} > ${CHOPPED_GENOME}

Then assign the midpoints of bins inside a gene, or to the nearest of any particular set of intergenic categories.

Here's a script I wrote which uses the output of closest-features to do that. This could be easily modified for any set of distance categories:

#!/usr/bin/env python

import sys
import numpy as np

'''
Input should be bin midpoints run through closest-elements, against some set of BED6+ (stranded) gene annotations
'''

def assign_to_category(dist):
  if dist < -50000:
    return '-50kb_to_-Inf'
  elif dist < -20000 and dist >= -50000:
    return '-20kb_to_-50kb'
  elif dist < -10000 and dist >= -20000:
    return '-10kb_to_-20kb'
  elif dist < -5000 and dist >= -10000:
    return '-5kb_to_-10kb'
  elif dist < -1000 and dist >= -5000:
    return '-1kb_to_-5kb'
  elif dist < 0 and dist >= -1000:
    return 'TSS_to_-1kb'
  elif dist == 0:
    return 'Gene'
  elif dist > 0 and dist <= 1000:
    return 'TES_to_+1kb'
  elif dist > 1000 and dist <= 5000:
    return '+1kb_to_+5kb'
  elif dist > 5000 and dist <= 10000:
    return '+5kb_to_+10kb'
  elif dist > 10000 and dist <= 20000:
    return '+10kb_to_+20kb'
  elif dist > 20000 and dist <= 50000:
    return '+20kb_to_+50kb'
  elif dist > 50000:
    return '+50kb_to_+Inf'

for in_line in sys.stdin:
  (bin_str, upstream_gene_str, upstream_gene_dist, downstream_gene_str, downstream_gene_dist) = in_line.rstrip().split('|')
  try:
    upstream_gene_dist = int(upstream_gene_dist)
  except ValueError:
    upstream_gene_dist = np.nan
  try:
    downstream_gene_dist = int(downstream_gene_dist)
  except ValueError:
    downstream_gene_dist = np.nan
  upstream_gene_elems = upstream_gene_str.split('\t')
  downstream_gene_elems = downstream_gene_str.split('\t')

  '''
  We have six possible scenarios, depending on the minimum distance and strandedness of closest elements
  '''

  '''
  1) If either distance value is NaN, the bin is assigned to the non-NaN distance element
  '''
  if upstream_gene_dist is np.nan and downstream_gene_dist is np.nan:
    raise SystemError("Bad NaN input on:\n[{}]\n".format(in_line.rstrip()))
  elif upstream_gene_dist is np.nan:
    gene_id = downstream_gene_elems[3]
    if downstream_gene_elems[5] == '+':
      # bin is closer to downstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      sys.stdout.write('{}\t{}\t\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue
  elif downstream_gene_dist is np.nan:
    gene_id = upstream_gene_elems[3]
    if upstream_gene_elems[5] == '+':
      # bin is closer to upstream TES
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to upstream TSS
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    continue

  '''
  2) If either distance value is zero, the bin is assigned to "Gene" to state overlap with a gene
  '''
  if upstream_gene_dist == 0:
    gene_id = upstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  if downstream_gene_dist == 0:
    gene_id = downstream_gene_elems[3]
    sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(0), gene_id))
    continue

  '''
  3) If the upstream strand is positive, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  4) If the upstream strand is negative, if the downstream strand is positive, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TSS
   -- the downstream distance is a signed (+) distance to the downstream gene's TSS

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '+':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TSS
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(downstream_gene_dist))), gene_id))
    continue

  '''
  5) If the upstream strand is positive, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '+' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TES
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(upstream_gene_dist)), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  6) If the upstream strand is negative, if the downstream strand is negative, then:

   -- the upstream distance is a signed (-) distance to the upstream gene's TES
   -- the downstream distance is a signed (+) distance to the downstream gene's TES

  '''
  if upstream_gene_elems[5] == '-' and downstream_gene_elems[5] == '-':
    if abs(upstream_gene_dist) < abs(downstream_gene_dist):
      # bin is closer to upstream TSS
      gene_id = upstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(-(abs(upstream_gene_dist))), gene_id))
    else:
      # bin is closer to downstream TES
      gene_id = downstream_gene_elems[3]
      sys.stdout.write('{}\t{}\t{}\n'.format(bin_str, assign_to_category(abs(downstream_gene_dist)), gene_id))
    continue

  '''
  Sanity check: If we get here, something went wrong
  '''
  raise SystemError("Bad logic on:\n[{}]\n".format(in_line.rstrip()))

To use it, e.g.:

$ MP_LEFT=500
$ MP_RIGHT=-499
$ MP_LEFT_UNDO=-500
$ MP_RIGHT_UNDO=499
$ CATEGORIZED_BINS=${PWD}/../data/categorized_bins.${ASSEMBLY}.bed
$ ./assign_chops.py < <(bedops --range ${MP_LEFT}:${MP_RIGHT} --everything ${CHOPPED_GENOME} | closest-features --dist - ${GENCODE_GENES}) | sort-bed - | bedops --range ${MP_LEFT_UNDO}:${MP_RIGHT_UNDO} --everything - > ${CATEGORIZED_BINS}

The output from this file annotates bins with the closest feature to a gene. It looks something like this:

$ more categorized_bins.hg38.bed
...
chr1  32000 33000 TES_to_+1kb RP11-34P13.3
chr1  33000 34000 TES_to_+1kb FAM138A
chr1  34000 35000 Gene  FAM138A
chr1  35000 36000 Gene  FAM138A
chr1  36000 37000 Gene  FAM138A
chr1  37000 38000 TSS_to_-1kb FAM138A
chr1  38000 39000 -1kb_to_-5kb  FAM138A
...

For example, the region chr1:33000-34000 is closest to the 1kb window downstream of the TES of FAM138A.

Once you have this file, you can quickly map it to your sorted regions-of-interest file:

$ bedmap --echo --echo-map --delim '\t' regions-of-interest.bed categorized_bins.hg38.bed > answer.bed

The file answer.bed will contain your regions of interest and the annotated bins which overlap them.


PCA on genotype matrix with multiple alleles

You can make a 'dummy variable' for each allele. That means that you don't have info per SNP, but for SNPs with more alleles, the allele is present (1) or not (0).


How to find common sequences among 6 multi-fasta files

If both the ID lines and the sequences are duplicated, you can do this using the FastaToTbl and TblToFasta scripts I have posted before and just normal Linux utilities:

FastaToTbl file1.fa file2.fa ... file6.fa |
    sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta > dupes.fa

To understand this, let's look at the output of each step. The FastaToTbl script
just converts fasta to a tab-delimited format where the ID is the first field and the sequence is the second, all on one line:

$ FastaToTbl file*fa | head -n1
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Then, we pass this through sort and uniq -dc which prints out the number of times a line was repeated (only for duplicated lines):

$ FastaToTbl file*fa | sort | uniq -dc
      4 some other sequence ALASOOPAKASKAPPPQQCDDREETAYSSLDAS
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Since we only want those that were present in all 6 files, we can use awk telling it to only print lines whose first field was 6:

$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6'
      6 Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

But then we don't want the 6 as part of the final output, so we need to remove it:

$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //'
Haladaptatus sp. R4 MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAESPGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTDAERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGRGTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTESALEEVADR

Finally, convert to fasta again:

$ FastaToTbl file*fa | sort | uniq -dc | awk '$1==6' | sed 's/^  *6 //' | TblToFasta
>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR

Alternatively, if your files are small enough to be loaded into memory, you could do this instead:

FastaToTbl file*fa |
    awk '{a[$0]++}END{for(s in a){if(a[s]==6){print s}}}' |
        TblToFasta > dupes.fa

Or you could just do the whole thing in a single awk operation:

awk '{
        if(/^>/){
            a[seq]++;
            seq=$0;
        }
        else{
            seq=seq"\n"$0
        }
      }
      END{
        a[seq]++;
        for(s in a){
            if(a[s]==6){
                print s
            }
        }
       }' file*.fa


Plot a circos plot to show the consistency between 2 samples

A circos plot is most likely not the most appropriate solution here. What I would suggest is a confusion matrix, of which you can find an example here:

For every variant in your vcf you'll add a number in this matrix. One sample is the columns, the other is the lines. If your variant is homozygous in both, then you add in that square +1 (the cell with 5845 in the example).

A perfect concordant sample pair will have only variants on the diagonal.

Here is some python code to get such a matrix. It uses cyvcf2 and pandas, and expects as input a vcf file with both samples.

from argparse import ArgumentParser
from cyvcf2 import VCF
import pandas as pd

def main():
    args = get_args()
    confusion_matrix(args.vcf)

def confusion_matrix(vcff):
    """
    First level of the dict is the "first" call, second level is the "second" sample
    0: hom_ref
    1: heterozygous
    2: unknown/nocall
    3: hom_alt
    """
    zygosities = {0: {0: 0, 1: 0, 2: 0, 3: 0},
                  1: {0: 0, 1: 0, 2: 0, 3: 0},
                  2: {0: 0, 1: 0, 2: 0, 3: 0},
                  3: {0: 0, 1: 0, 2: 0, 3: 0},
                  }
    for v in VCF(vcff):
        zygosities[v.gt_types[0]][v.gt_types[1]] += 1
    zygs = [2, 0, 1, 3]
    df = pd.DataFrame(index=zygs, columns=zygs)
    for tr in zygs:
        for te in zygs:
            df.loc[tr, te] = zygosities[tr][te]
    df.columns = ['nocall', 'hom_ref', 'het', 'hom_alt']
    df.index = ['nocall', 'hom_ref', 'het', 'hom_alt']
    print(df)

def get_args():
    parser = ArgumentParser(description="Create confusion matrix of SNV calls")
    parser.add_argument("vcf", help="vcf containing two samples")
    return parser.parse_args()

if __name__ == '__main__':
    main()


Off-by-one error when mapping Ensembl IDs to gene symbols

There’s no 1:1 mapping between the IDs of different database schemas. As a consequence, some Ensembl gene IDs map to multiple MGI symbols, or to none (and vice versa).

Therefore, you can’t assume that the number of result rows will be equal to the number of query IDs.


Importing a WGCNA co-expression network into a networkX graph in Python

You can use exportNetworkToCytoscape or exportNetworkToVisANT to generate and egde list file with the two nodes at the ends and weight for each edge. Hopefully there's a way to import that into networkX.


How to change sequence format in my alignment file?

This is a comment rather than answer - but not enough rep

I don't think you'll find any tool to do this exactly as you want  - how is it going to decide what's a dot and what's a letter (or a space,N or whatever) if you have multiple ambiguities or even an even number of alignments split 50:50?

Wouldn't be hard to code something which prints the whole column if it consists of non-identical nucleotides, else dots. Would be easy in R, could probably even do it in Excel without too much effort and no need to use VBA.


How I can test my hypothesises computationally

This is what packages like FateID and Monocle are for, namely taking single-cell RNAseq data and inferring differentiation trajectories from it. Don't try to reinvent the wheel on this, there are a number of packages out there to do this sort of thing and they're going to have pretty complicated methods that you're really not going to want to reinvent unless you absolutely have to.


Most scalable pseudotime ordering algorithm

As it turns out, pseudotime algorithms are exceedingly difficult to parallelize. Pseudotime does not scale linearly to a subset of data compared to the full dataset, and as such a course-grained approach to parallelizing pseudotime calculations on random bins of the data, followed by bin merging, does not work. These values will never converge on the actual full-dataset solution.

I put in a pretty good effort trying to parallelize slingshot, but with no success. I realize my code may be of limited utility without more explanation, but I'm not putting more time into this. FWIW:

This code:

Finds computationally optimal bin size for slingshot based on a number of genes and features
Randomly splits the SingleCellExperiment object into that number of bins and runs slingshot on each bin
Orients each bin so that a cell with a pseudotime of 1 in the first bin is similar to a pseudotime of 1 in every bin. Pseudotime can get flipped depending on the sampled subset.
Merges all bins into a single matrix so each cell has a assigned pseudotime based on randomly sampled bins.
Repeats steps #2 to #5 with a different random seed for binning every time.
Measures convergence towards actual solution at three pseudotime values, where convergence is defined as correlation of cell indices along the trajectory.

I abandoned the project after observing that more iterations to not change the convergence towards the actual solution. While pseudotime values close to 0 recapitulated the actual solution quite well, pseudotime values towards 100 failed to correlate at all with the actual solution. There was no convergence with more iterations.

Maybe this will be helpful to somebody, which is why I'm posting.

# Parallel Iterative Sampling with Slingshot
# pissshot is a course-grained estimator of the actual slingshot solution, and iteratively converges towards the actual solution
# pissshot randomly partitions cells in a large database into bins, bin size is determined by finding the computationally optimal number of cells
# pissshot is run independently on each bin, pseudotime vectors are aligned across all bins, and then all bins are merged
# pissshot repeats this process and the average of pseudotime assignments for each cell is taken at the end of each iteration.
# The convergence of the blingshot model at the end of each iteration is calculated by running slingshot on a small sample of adjacent cells at a random pseudotime interval and measuring correlation with the predicted model
# When satisfactory convergence is reached, blingshot returns a singlecellobject, just as slingshot would

pissshot <- function(sce, num_iterations = 20){
    cat("Running blingshot on", dim(sce)[2],"cells and",dim(sce)[1],"features\n")
    cat("  Step 1/5: Finding computationally optimal increment sizes\n   ")
#    increment_size <- FindOptimalIncrementSize(sce)
    increment_size <- 688
    num_increments <- ceiling(dim(sce)[2]/increment_size)-1
    cat("\n  ...Increment size of",increment_size,"cells is computationally optimal\n")
    cat("  Step 2/5: Course-grained pseudotime assignment across",num_increments,"increments")
    ps <- list()
    convergenceArr <- c()
    for(iteration in seq(from=1,to=num_iterations,by=1)){
        cat("\n   ... iteration",iteration,"...\n")
        cat("         randomly assigning cells to",num_increments,"bins\n         0%.")
        bins <- PartitionBins(sce, increment_size = increment_size, seed = 7*iteration)
        cat("100%\n         running slingshot on each bin\n         0%")
        cell_pseudotimes <- c()
        for(b in seq(from=1,to=length(bins),by=1)){
            bin_with_pseudotime <- suppressMessages(slingshot(bins[[b]]))
            # Find right cell and left cell in pseudotime array (cell with highest/lowest pseudotime values)
            pcells <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1)
            rownames(pcells) <- colnames(bin_with_pseudotime)
            pcells <- data.frame(lapply(pcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(pcells))
            pcells.sorted <- pcells[order(-pcells[,1]), , drop = FALSE]
            right_cell_ID <- rownames(pcells.sorted)[dim(pcells.sorted)[1]]
            left_cell_ID <- rownames(pcells.sorted)[1]
            right_cell_pos <- match(right_cell_ID,colnames(bins[[b]]))
            left_cell_pos <- match(left_cell_ID,colnames(bins[[b]]))
            right_cell_logcounts <- assays(bins[[b]])$logcounts[,right_cell_pos]
            left_cell_logcounts <- assays(bins[[b]])$logcounts[,left_cell_pos]
            if(b==1&&iteration==1){
                ref_right_cell_logcounts <- right_cell_logcounts
                ref_left_cell_logcounts <- left_cell_logcounts
            }
            # Figure out whether this pseudotime in this bin has flipped relative to reference (first bin), and if so, flip pseudotime values
            cor_opposite_ends <- mean(cor(right_cell_logcounts,ref_left_cell_logcounts),cor(left_cell_logcounts,ref_right_cell_logcounts))
            cor_same_ends <- mean(cor(right_cell_logcounts,ref_right_cell_logcounts),cor(left_cell_logcounts,ref_left_cell_logcounts))
            flipped <- FALSE
            if(cor_opposite_ends > cor_same_ends){
                flipped <- TRUE
            }
            if(flipped==TRUE){
                # find max pseudotime in colData(bin_with_pseudotime)$slingPseudotime
                # Recalculate pseudotime in bin_with_pseudotime as max-value
                pseudotimes <- colData(bin_with_pseudotime)$slingPseudotime_1
                colData(bin_with_pseudotime)$slingPseudotime_1 <- max(pseudotimes)-pseudotimes
            }
            cells_with_pseudotime <- as.data.frame(colData(bin_with_pseudotime)$slingPseudotime_1,colnames(bin_with_pseudotime))
            cell_pseudotimes <- rbind(cell_pseudotimes, cells_with_pseudotime)
            cat(".")
        }
        cat("100%\n")

        ps[[iteration]] <- cell_pseudotimes[order(row.names(cell_pseudotimes)), , drop = FALSE]

        pss <- data.frame(ps)
        if(iteration>1){
            convergence10 <- measureConvergence(pss, sce = sce, ptime = 10, increment_size = increment_size)
            convergence50 <- measureConvergence(pss, sce = sce, ptime = 50, increment_size = increment_size)
            convergence99 <- measureConvergence(pss, sce = sce, ptime = 99, increment_size = increment_size)
            cat("   convergence at 10:",convergence10,"... at 50:",convergence50,"...at 100:",convergence99,"\n")
            convergenceArr <- rbind(convergenceArr,c(iteration,convergence10,convergence50,convergence99))
        }
        # Iterate until the model converges towards the actual slingshot solution, within an indicated fraction (i.e. convergence = 0.02)
        # The convergence measure compares how similar the model is to the actual slingshot solution, i.e. a convergence of 0.02 means that cell pseudotime values are within 2% of the actual solution
        # Assess this by selecting a pseudotime at random, analyzing the following increment_size cells, and comparing accuracy of pseudotime measurement across that trajectory with the composite of the course-grained analysis

        # need to study how convergence changes over pseudotime
    }

    # remove cells with really high standard deviations

    return(pss)
#    return(convergenceArr)
}

measureConvergence <- function(mat, sce = sce, ptime = 45, increment_size = 500){
    # this line is throwing an error
    mat[,"avg"] <- apply(mat[,1:dim(mat)[2]],1,mean)
    mat <- mat[order(mat[,"avg"]), , drop = FALSE]
    # select a random pseudotime increment
    mat <- mat[mat[, "avg"] > ptime,]
    mat <- mat[1:increment_size,]
    cell_IDs <- rownames(mat)
    # pull out this list of cell_IDs from the sce object and run slingshot
    sce2 <- suppressMessages(slingshot(sce[,cell_IDs]))
    cpcells <- as.data.frame(colData(sce2)$slingPseudotime_1)
    rownames(cpcells) <- colnames(sce2)
    cpcells <- data.frame(lapply(cpcells, function(x) as.numeric(as.character(x))), check.names = F, row.names = rownames(cpcells))
    # sort cpcells by increasing pseudotime
    cpcells <- cpcells[order(cpcells[,1]), , drop = FALSE]
    ordered_cells <- cbind(rownames(mat),rownames(cpcells))
    # get the index of each cell in the dataframe, run a correlation on how well the indices line up
    cell_indices <- c()
    for(cell in cell_IDs){
        new_indices <- c(which(ordered_cells[,1] == cell), which(ordered_cells[,2] == cell))
        cell_indices <- rbind(cell_indices,new_indices)
    }
    convergence_measure <- cor(cell_indices[,1],cell_indices[,2])
    plot(cell_indices[,1],cell_indices[,2])
    return(convergence_measure)
}

PartitionBins <- function(sce, increment_size = 500, seed = 123){
    bins <- c()
    while(dim(sce)[2]>=increment_size){
        if(dim(sce)[2]<increment_size*2){
            bins <- c(bins, sce)
            sce <- sce[,1]
        }
        else{
            set.seed(seed)
            train_ind <- sample(seq_len(ncol(sce)), size = increment_size)
            bins <- c(bins, sce[,train_ind])
            sce <- sce[,-train_ind]
            cat(".")
        }
    }
    return(bins)
}

FindOptimalIncrementSize <- function(sce){
    # Find the computationally optimal increment_size for grainy alignments using 5 increment points (100, 500, 1000, 2500)
    res <- c(100,250,500,750,1000,1500,2500,5000)
    runtimes <- c()
    for(i in res){
        runtimes <- c(runtimes, suppressMessages(system.time(slingshot(sce[,0:i])))["elapsed"])
    }
    mat <- as.data.frame(t(rbind(res,runtimes)))
    c <- summary(lm(mat$runtimes ~ poly(mat$res, 2, raw=TRUE)))$coefficients[,"Estimate"]
    predicted_runtimes <- c()
    for(i in seq(from=100, to=5000, by=1)){
        predicted_runtime <- (c[1] + c[2]*i + c[3]*i^2)*(dim(sce)[2]/i)
        predicted_runtimes <- rbind(predicted_runtimes,c(i,predicted_runtime))
        cat(".")
    }
    increment_size <- predicted_runtimes[which.min(predicted_runtimes[,2]),1]
    return(increment_size)
}


Store and search on variant and annotation data

The question you linked is about storing genotypes. Your question is about keeping site information, which is different and much easier.

As to solutions, you can parse the information you want and store them in a standard SQL database. You can use UCSC binning for region queries. Modern SQL databases all provide R-tree indices, which can be adapted for region query and is probably better these days.


Expand granges object different amounts upstream vs. downstream

There doesn't appear to be a built in function for that, so here's a modified version of promoters():

expandRange = function(x, upstream=2000, downstream=1000) {
  strand_is_minus = strand(x) == "-"
  on_plus = which(!strand_is_minus)
  on_minus = which(strand_is_minus)
  start(x)[on_plus] = start(x)[on_plus] - upstream
  start(x)[on_minus] = start(x)[on_minus] - downstream
  end(x)[on_plus] = end(x)[on_plus] + downstream
  end(x)[on_minus] = end(x)[on_minus] + upstream
  x
}

You can then call expandedGenes = expandRange(genes, 5000, 1000).


Retrieving all genes for a gene ontology term

It is filtering by your GO terms, it is only giving you genes associated with those GO terms. However the GO terms in the output are not linked to the GO terms in your input, so it is giving you all the GO terms linked to the genes which are associated with the GO terms you filtered by.
Essentially BioMart is not very clever. It has a table of all the GO terms linked to the genes, it filters the genes to only include those linked to the GO terms you've given it. Then it takes that list of genes and queries that same table completely separately, getting the GO terms for each gene.


Yes, your understanding is largely correct. This originates from the situation that for each detector on a nanopore array there are 4 pores. I'll explain mux scans and groups, but this is outdated information as now another system is used.

So a MinION FC has 2048 pores with 512 sensors. Not all of these pores will be equally suitable for sequencing. At the beginning of the run originally a mux scan was performed, to test for each sensor which the best pore would be (if any) out of the 4 available. The best pores will be used for sequencing first, the rest will be turned off. That is group 1. After 8 hours the expectation is that many of these pores have degraded/accumulated damage or are not functional for other reasons. Then it is time to swith to group 2. So you change to another set of pores, which in these cases often lead to a temporary increase in yield. The default script used 8 hours, although some used adaptive software to select new pores when the first group started going bad. Also note that during the run no new mux scan was performed: the software did not check if the selection it made in the beginning still made sense. Pores in group 1 got turned off after 8 hours and switched to group 2, without taking their current 'health' into account. Cumulative yield plots get a characteristic bumpy pattern from this. Each bump is less than the previous one, eventually flatting off.

Now that's history, as the software now uses dynamic muxing and will continuously check for the best pore per channel - and as such maximize yield by having as many healthy pores as possible available for sequencing. Cumulative yield plots are now smoother.


Minimal working example for phylogenetic tree construction in R using an alignment of DNA strings and a vector of names?

A key difference from many examples provided is that they assume that the data provided is in the form of sequences rather than an already aligned sequence alignment.
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

#read in the data from a format with the aligned sequences (eg msf)
msf.res <- read.alignment(file = "geneticAlignment.msf", format = "msf")
#dist.alignment and not dna.dist commonly used in tutorials
msf.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
#perform neighbor joining upon the distance set
msf.res.dist.alignment.nj = nj(msf.res.dist.alignment)
#then plot the joined sequences
plot(msf.res.dist.alignment.nj, main = "Neighbor Joining")


I don't know if you have followed the discussions from the last years, but FPKM and RPKM are obsolete. Statisticians found that FPKM and RPKM are not the right values that can be used for proper statistics. It is a wrong way of normalizing RNA-seq data. In stead they suggest to use raw read counts, which can be used in tools like edgeR or DESeq2. Also, CPM would be better than FPKM or RPKM, for example you can use (log2) CPM values in limma with trend.

But to answer your questions.

You can't convert, you need to recalculate from raw data.
No not acceptable. To compare studies you need to correct for batch, which you can do with statistical programs such as edgeR and DESeq2.

The best advice I can give is to start with raw data, also from the already published data that you would like to compare your own experiment with. And continue with edgeR or DESeq2 or the like.


It is actually very hard to answer this using only the information you have provided. It'd be necessary to actually see your test.nf file and read the processes. The execution seems right and I'd guess that the processB may be not emiting any output as shown in https://www.nextflow.io/docs/latest/process.html#outputs.
It is also possible that when using the -resume together with a specific run-name you may be recovering the state of that specific run as shown here https://www.nextflow.io/blog/2019/demystifying-nextflow-resume.html. You can try to execute it without the run-name or without the -resume flag so you can check whether the problem is with the emition and channelling of outputs between processes, or with the -resume flag.


No refseq transcripts in 2nd half of Y chromosome for HG38

Nevermind, Answered my own question. It's a repetitive region that isn't very well understood

  The human Y chromosome, because it is enriched in repetitive DNA, has been very intractable to genetic and molecular analyses. There is no previous evidence for developmental stage- and testis-specific transcription from the male-specific region of the Y (MSY). Here, we present evidence for the first time for a developmental stage- and testis-specific transcription from MSY distal heterochromatic block. We isolated two novel RNAs, which localize to Yq12 in multiple copies, show testis-specific expression, and lack active X-homologs. Experimental evidence shows that one of the above Yq12 noncoding RNAs (ncRNAs) trans-splices with CDC2L2 mRNA from chromosome 1p36.3 locus to generate a testis-specific chimeric β sv13 isoform. This 67-nt 5′UTR provided by the Yq12 transcript contains within it a Y box protein-binding CCAAT motif, indicating translational regulation of the β sv13 isoform in testis. This is also the first report of trans-splicing between a Y chromosomal and an autosomal transcript.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1832090/


How to concatenate two fasta sequences by py

First of all, those are not fasta sequences. Fasta is a specific format, and requires a > followed by a sequence identifier on one line and the one or more lines of sequence.
Assuming you really have two files in the format you show, you can do:
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]; next} for(b in bbone){print b$2}}' file1 file2
GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

The gsub(/ /,""); just removes the extra space you have in file1 after each :. If you also want to have the sequence names as you show, use:
$ awk -F: '{ gsub(/ /,""); if(NR==FNR){bbone[$2]="Backbone"NR; next} for(b in bbone){print bbone[b]"+"$1":", b$2}}' file1 file2
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCG
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG

If you really need to do this in python for some reason, please edit your question and show us your python script so we have something to work on and expand. I really don't want to write a whole python program for something so trivial.


As far as I can tell, you use the -lib option for species where there is no existing RepeatMaster database. Therefore, all it does is attempt to locate versions of the repetitive elements from the -lib file in the target genome file. However, since your lib file doesn't have any information on what type of element each of your hits is, it looks like RepeatMasker is treating all of them as "Unclassified".

Note the following from the repeatmasker.help file you should have downloaded along with RepeatMasker:

  The recommended format for IDs in a custom library is:

  >repeatname#class/subclass

  or simply

  >repeatname#class

  In this format, the data will be processed (overlapping repeats are
  merged etc), alternative output (.ace or .gff) can be created and an
  overview .tbl file will be created. Classes that will be displayed in
  the .tbl file are 'SINE', 'LINE', 'LTR', 'DNA', 'Satellite', anything
  with 'RNA' in it, 'Simple_repeat', and 'Other' or 'Unknown' (the
  latter defaults when class is missing).

This makes sense. RepeatMasker cannot know if an arbitrary sequence you give it is a SINE or a LINE or anything else. It will therefore rely on the categories present in the fasta headers. Since your file has no categories, it identifies everything in your lib file as "Unknown", which is presumably what the "Unclassified" section in your output refers to:

Unclassified:    560941     66374119 bp   31.69 %

Total interspersed repeats: 66374119 bp   31.69 %


Error in creating a volcano plot in MATLAB

I'm not sure if nowadays matlab is still a good option for informatics. I would rather go for a R script to make a volcano plot. If you don't want to have fancy code but functional one, you can even make the foldchange calculations and log2 transformation in excel and just plot the log2 fold change in x axis and the -log10 transform of the p/q values on the y axis.
If you need it, I could provide one of my own scripts.

Update: here is the actual script:

#load data into the 'res' variable
res <- read.table("genes.txt", header=TRUE)
head(res)

# Make a basic volcano plot with grey dots
with(res, plot(log2FoldChange, -log10(qvalue), pch=19, col="#4A4E4D", xlab="Log2 Fold Change", xlim=c(-8,8)))

# Add colored points: "indianred" if padj<0.0001, orange of log2FC>2, green if both)
with(subset(res, qvalue<.0001 ), points(log2FoldChange, -log10(qvalue), pch=19, col="#17B090"))
with(subset(res, abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#F89538"))
with(subset(res, qvalue<.0001 & abs(log2FoldChange)>2), points(log2FoldChange, -log10(qvalue), pch=19, col="#FA6367"))


Create GFF3 feature exporter. The information is below:

Your answer looks like a good start. You'll probably want to use sys.argv, argparse, or another argument parsing tool to read in options from the command line, then compare them to the fields (e.g. the "type" argument gets compared with the type column).
Given that attributes can be arranged in an arbitrary order, you'll probably need to read all of them from each line and store in a dict in order to determine whether the target is present. A regular expression search is another way to do that.
There are a few libraries available (e.g. from biopython) for reading and subsetting FASTA records.


Basic question about finding mRNA sequence in transcriptome

Basically you want to know the VHL Gene expression in that 10x dataset. They have provided the barcode, matrix and feature files on GEO GSE131685 as Supplementary file. You could use those files and get the gene expression matrix as described on 10x website.
library(Matrix)

matrix_dir = "/opt/sample345/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")

mat <- readMM(file = matrix.path)

feature.names = read.delim(features.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

barcode.names = read.delim(barcode.path,
                           header = FALSE,
                           stringsAsFactors = FALSE)

colnames(mat) = barcode.names$V1
rownames(mat) = feature.names$V1

This mat can be further loaded into Seurat for any downstream analysis. All this can be done on a laptop.
If you want to re-quantify the gene expression for whatever the reason, you need to re-run the cellRanger pipeline using GTF of your interest (it needs some preparation), which may not work on your laptop. It requires large memory and long running times, ideally to submit as job on cluster.


Perhaps,removing introns (if it is a higher eukaryote)

extractfeat tembl:x65921 -type CDS -join

It is difficult to know without knowing the gene.

BTW NCBI's Entrez "efetch" will do essentially the same thing, but its more complicated


Alternative to enrichR for enrichment analysis?

A couple of years ago I wrote a blog post about GSEA methods in Bioconductor, in addition to these I would recommend clusterProfiler or ReactomePA.  In most of them you need to provide the gene set/pathway database.
If you want to do gene ontology enrichment then I recomment the topGO package.
In most, if not all, of these methods you just need an ordered list and a gene set source.


How to reproduce this specific sequence logo from 10 000 human mRNAs?

If you look at the Appendix from https://www.ncbi.nlm.nih.gov/pmc/articles/PMC306349/, there is a table of sequences that were likely used to construct the sequence logo, at least partially, from the 5' end to the stereotyped ATG subsequence.

It's not in a nice format amenable to use directly with weblogo (like FASTA), but with some programming/elbow grease it might be rendered into that form, which you can then use.

Generally, you might investigate the supplementary materials of these and related papers, to see if you might something directly usable.


Interpretion of my coronavirus 2019-nCov, Wuhan, China BLAST tree?

The  evolutionary related group (clade) of betacoronaviruses you have identified share an amino acid homology of 85% and include SARS. I know this from the underlying tree published on BioRxiv of a broader group of betacoronaviruses, i.e. your data is a defined subset of the betacoronaviruses which all share a unique, single common ancester.
Lets call this group the SARS clade.
You have performed a nucleotide blast and asked NCBI to produce a NJ tree using 2019-nCov as the reference. I can tell this by the colour coding and the scale bar at the bottom right hand corner shows the genetic distance involved is reasonably larger than 15% divergence by amino acid data. The scale bar represents the number of mutations per nucleotide.
In summary, your tree is essentially a subset of the broader genetic diversity of the one published, but there is a rooting issue.
In your tree, the majority of sequences are from the 2002 SARS epidemic and the virtually zero genetic distance between them is simply because it is a rapidly transmitted outbreak. I didn't realize that SARS has two independent origins both initially from bats. This is quite scary.

2019-nCov is an outgroup within the 'SARS clade' hence it appears at the other side of the tree, i.e. they share a more distant common ancestor
However 2019-nCov is not the most distant common ancester within this subset of the betacoronaviruses, i.e. SARS clade, this belongs to the two bat viruses ZC45 and ZX21. The program has made a likely rooting error (below).
Again the blast omitted the majority of the betacoronaviruses, for example MERS

Rooting issue the reason I suggest there is a rooting error is because the BioRxiv tree using a broad sample of the betacoronaviruses placed the bat strains ZC45 and ZX21 as outgroups to the SARS clade and 2019-nCOV was immediately within that, so in this definition 2019-nCOV is an 'ingroup' within the SARS clade, whereas your tree it is an 'outgroup'. It not a huge issue, but the location of the "root" is determined by the common ancester above it (sequences with <85% homology) and in your tree those have been omitted.
Generally, I like the analysis otherwise and gives an insight into SARS that I hadn't previosly been aware of.


Machine learning using protein-sequences

There are two types of features which can be extracted from protein sequences.  It is not possible to known which, if any, of these features would be useful in classification tasks.  It may not be possible to build a classifier at all, or it may be very straight forward.  To know this, a feature selection technique must be used, such as e.g. forward feature selection.

Sequence based features

Sequence based features are derived from sequence-only information.

The most common used features are the PseAAC (pseudo amino acid composition) style features, including:

Average Position of each amino acid - 20 features
Average Distance between each amino acid - 20 features
Motifs - 20 features, 20*20 features, and/or 20*20*20 features (depending on motif length)
OAAC (Overall Amino Acid Composition) - 20 features
SAAC (Split Amino Acid Composition) - 60+ features
AAIndex Composition - 1 to 567+ features
PCC (Physicochemical Composition) - ~13 features
PSSM encoded as 1, 20, 200, 400 or 600 features (i.e. Psi-Blast output using NR or swissprot)
Average predicted residue burial/exposure level
Average predicted sequence entropy

Structure based features (predicted)

Using the sequence, you can either predict the 3d structure using homology, ab initio, or directly use a homologous structure.  You can then extract predicted 3d features, such as:

3d size (i.e. max of x,y,z)
Calculated ASA/SASA
3d contacts (type, frequency, etc)
any other local/global environmental information
cath/scop classification
dssp/stride composition


STAR Indexing Diference for Small and Large Genome File Output

You didn't get extra information with the small genome, the index for the larger one just wasn't made. My guess is that you ran out of RAM. For mammalian sized genomes ensure you have at least 30GB of RAM, likely more like 60GB. Also, make sure to look for error messages.


While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError'

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message.


How to count reads in bam per bed interval with bedtools

The order of -a and -b switched at some point. You want:

bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt

For reference, this is the end of the help output in version 2.25:

Default Output:
     After each entry in A, reports:
       1) The number of features in B that overlapped the A interval.
       2) The number of bases in A that had non-zero coverage.
       3) The length of the entry in A.
       4) The fraction of bases in A that had non-zero coverage.

And this is the equivalent output from version 2.19:

Default Output:
     After each entry in B, reports:
       1) The number of features in A that overlapped the B interval.
       2) The number of bases in B that had non-zero coverage.
       3) The length of the entry in B.
       4) The fraction of bases in B that had non-zero coverage.

Update: The change in behavior happened in version 2.24:

  We have changed the behavior of the coverage tool such that it is consistent with the other tools. Specifically, coverage is now
  computed for the intervals in the A file based on the overlaps with the B file, rather than vice versa.


What is "S13"?

S# is most likely the sample number. If so, you do not want to combine these.

What is "L002_l1"?

L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).

What is "L002_R1"?

R1 is one end of paired-end sequencing, R2 is the other end. These filepairs must have the same number of lines. Use an aligner that utilizes both.


How is the t-statistic value calculated in GEO2R or Bioconductor?

That is not trivial. GEO2R uses the Bayesian linear model-based framework limma https://bioconductor.org/packages/release/bioc/html/limma.html for the analysis. I doubt you can easily and meaningfully implement the relevant code in Python without a major effort and notable background knowledge. The user guide https://bioconductor.org/packages/release/bioc/vignettes/limma/inst/doc/usersguide.pdf has background information plus links to the relevant papers. Don't try to analyze things by hand, use expert software. If it is not significant, well, then it is not significant. A 2vs2 comparison is strongly underpowered, therefore lack of DEGs are not unexpected, especially for microarrays.


What is "aligned sequences" and "consensus sequence" in the context of sequence logo? How to compute these?

Panel (a) shows a sequence gap-free alignment. Each row corresponds to a contiguous 15 base pair sequence of DNA (e.g. row 1 could be a human sequence, row 2 could be the equivalent mouse sequence etc.). Each column corresponds to a particular position/residue (e.g. column 1 is the 1st position of the alignment etc.).

Panel (b) shows the sequence logo. As noted in the comments, this does not correspond to the alignment in (a). Instead, here is the corresponding logo:

From this you can intuitively understand that at the 1st sequence position you have only As and Ts and their frequency corresponds to the y-axis. For the 3rd position, you have only As, so only A is represented here in the logo etc.

The consensus sequence is the most frequent letter at each position in the alignment, i.e. [AT]TAGC etc. where [AT] means A or T since they are jointly high in frequency (5 each).


RNAseq DE comparison with samples of different read length

You don't need to do anything, this is taken care of by the linear model - if you don't have any confounding effects.
If you do have confounding effects, e.g. all control samples on HiSeq and all treated on NovaSeq there's not a lot you can do anyway, other than be honest when you present your results.
Taken care of by the model means adding an extra 2 level factor for sequencer by sample, you need to do this anyway to take care of the batch effect between sequencer runs.


Targeting Illumina reads with >70X depth-of-coverage is a good place to start.  Once a robust database is established lower coverage isolates can be used effectively too.
More specific to individual SNPs... filtering on VCF QUAL value, Map Quality and/or AC value have shown to be useful.
Core genome will likely be around 1,000 SNPs, but it will depend on the lineage the reference belongs to and the aligned sample lineage.  If both are from the same lineage you're likely to be less than 700 SNPs.  When from different lineages SNP counts will obviously be higher in addition to the alignment error seen.
Pertaining to filtering individual SNPs...
vSNP is a tool developed for SNP analysis.  It can be applied to different specie types but was originally developed for use with TB complex isolates.  Code detail is here.
vSNP is a 2-step process.  First, generating VCF files.  Second, using those VCF files to output SNP tables and trees.  The goal of the script is to achieve a high resolution genotyping analysis that can be validated with minimal effort.  Some steps are included to help achieve this such as reference selection, automatic grouping and position filtering.
Pseudocode:
For set of VCF files created from the same reference...
def find_positions()
    AC=2
    QUAL > 300

For each found position...
def get_snps()
    AC=2 (solid call)
    QUAL > 50
        THEN --> SNP
    AC=1 (mix call)
    QUAL > 150
        THEN --> AMBIGIOUS
    AC=2
    QUAL >= 50
       THEN --> REF
    QUAL < 50
       THEN --> N
    ALT == "None"
       THEN --> -

First find all possible positions that could be informative in all samples/VCF files: find_positions(), then define each position to be included in the alignment: get_snps().
Then it can be convenient to filter consistently poor alignment regions such as PPE/PGRS.


Are soft-clipped bases used for variant calling in samtools + bcftools?

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment. You'll note that the soft-clipped region isn't used in the depth/pileup (both tools use the same underlying code, so it doesn't matter which you use). If you're curious, samtools ignores soft-clipped bases because it's based on making a per-base stack of alignments covering each position. In the BAM format, alignments are sorted and assigned to bins according to their start/end positions, which won't include soft-clipping. Consequently, when samtools is making the pileup it won't even see the alignments that would overlap a given base if soft-clipped bases were included.

This then sort of begs the question of what GATK's HaplotypeCaller is doing differently. There, regions in the genome are essentially assembled in a small de Bruijn graph, which allows for soft-clipped bases around indels to then be resolved, given that the graph would start/end a little-way on past each side of indels. This is also why you don't need to do indel realignment with the HaplotypeCaller (this was needed in the old UnifiedGenotyper).

Edit: For more details regarding the HaplotypeCaller, see this nice page on GATK's website, which goes into much more detail than I did here.


Can I export the SNPs between whole genome alignments using the command line?

It works for me with the download from download page. The download archive contains Mauve.jar so the command line in the answer you link to should work.

However, you might need to include all the jar files (not just Mauve.jar) in the classpath for this to command work, though. For instance, the following works for me:

MAUVE_DIR=‹the-dir-containing-Mauve.jar›
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
java org.gel.mauve.analysis.SnpExporter -f aln.xmfa -o aln.snps

… there’s probably a smarter way to make this work but I don’t use Java enough to know how.


Why do BAM files created by different tools have different file sizes?

It's worth bearing in mind that when outputting compressed BAM, as most tools do by default, they may well be using different levels of compression and/or different libraries, or versions of said, libraries for doing (de)compression which will result in different file sizes.  Additionally coordinate sorted BAM will compress more than unsorted BAM. The current version of Picard uses HTSJDK which in turn uses java.util.zip.Deflater/Inflater, current versions of samtools should be using the HTSlib library which in turn depends on the standard zlib library. You can see the effect of different implementations of zlib have on file size and execution time in benchmarking done by the samtools team.

However, in your case the best way to see if there is any difference between the BAM files is to rule out the effect of different levels of compression or libraries used for compression and save both BAM files as uncompressed.  Both samtools and Picard have options for disabling or changing the levels of compression, since the BAM compression standard BGZF is implemented on top of the gzip format it has inherited the ability, just like with gzip, to change the compression level from 0 to 9.

samtools view -bu will allow you to produce uncompressed BAM output (which is also handy for piping into other programs as it saves time wasted compressing decompressing what is essentially a stream).  Also note that samtools sort has a -l INT setting where INT can be set between 0 (compression off, as with -u) 1 (for fastest compression, but increased file size) or -9 (for maximal compression, with increased run time).  Some of the effects of increased runtime for higher compression settings might be ameliorated using the -@ argument which allows you to set the number of extra threads used for BAM compression, by default samtools won't use any.

Picard tools has a general setting COMPRESSION_LEVEL which is applicable to most of its tools setting this to 0, COMPRESSION_LEVEL=0 should disable compression.

So re-running whatever Picard tool you used in the first instance with COMPRESSION_LEVEL=0 will then enable you to check that the file is not further modified by samtools view -bu. The assumption here being that if both files have exactly the same content they should be the same size uncompressed, of course if they have trivial differences with regard to white space formatting things may still differ.


What are the algorithms for biometric identification in large database?

I think that a very nice example of what you are talking about saying the word bigger, is probably a microarray transciptional analysis of RNA extracted from a biological sample. Simplifying, in these analyis, each observation is analised on, up to 10000 genes expression levels of different mRNA coding genes. Thus, from these type of molecular biology experiments, we obtain matrices with few observations (d, the samples) respect to the dimensions observed (k, the target gene expression levels). Well, d << k. In this situation, in order to perform a classification-like task, and to reduce the problem of the overfitting related with dimensionality maledition, the ways that we might follow are several:

1) PCA: it is a good tool to visualize multidimensional datasets reducing dimensionality and noise level in the consolidated and standardized data. Problem: it is an unsupervised machine learning algorithm. Alone, it haven't the capacity of class discrimination.

2) PCA-LDA is a good combination, where PCA eliminates  hyperdimensionality and background in our data sheets, and than LDA (linear discriminant analysis) finds the directions onto the new subspace of characteristics, that are able to maximize variances between classes and decrease the variances within these. For this reason LDA is known as a kind of supervised machine learning technique. Furthermore, PCA and LDA are usually used together because LDA applied on high dimensional datasets could lead us toward overfitting problems. After the training of PCA-LDA pipeline, being LDA a classificator object trained and shaped on the principal component subspace, you can use this to examine new unknown samples or observation and predict its class by similarity. e.g. recently, Martin et al., on a Nature article talk about PCA-LDA approach to identify and than to predict cancer cells from those are normal, starting from ranman spectra of tissues.

3) PCA-LDA on transpose matrix: to use as better as possible the huge dimensionality, some authors used to transpose the initial matrix in order to force positioning of genes as rows and samples as columns. Then, the subsequent LDA is performed not on PC scores, but on loading factors of the eigenfactor matrix obtaining thus the list of the most important discriminant gene expression levels to segregate normal samples to cancer ones.


Genomics: Downloading genomic data from the NCBI Genbank database

For downloading data from commandline use the E-utilities of NCBI.


Calculating read average length in a Fastq file with bioawk/awk

This can also be done with regular awk.

awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>

The NR%4==2 count the second line out of every block of 4.  length is a built-in that defaults to the length of the line, same as length($0).  In this case you can inject you custom printing to the END{} block but countread and countstored will always be the same since the averaging is done on the fly.


"Transforming" Fasta files to GTF or GFF3 files is a common request, as are similar tasks such as "converting" BAM files to VCF files. But these are underdefined questions, and "transform" or "convert" isn't really the right way to think about the task.

As far as your particular request is concerned, Fasta files and GTF files contain fundamentally different types data.

Fasta files contain sequences, which could be chromosomes, assembled contigs, polypeptides, transcripts, or other biological sequences.
GTF (and GFF3 and BED) files contain intervals, indicating that some feature of interest occurs at such-and-such position on such-and-such sequence. A hypothetical example: a transcription factor binding site occurs on chromosome 2 from position 40011 to 40017.

So you are correct: it's not really helpful to think of your task in terms of "transforming" Fasta to GTF. The process of creating a GTF/GFF3/BED file is often called annotation. There are many different methods to annotate many different types of genomic features: transposable elements, segmental duplications, protein-coding genes, non-coding RNAs, and so on.

As a result, it's meaningless to say "I need a GTF file" without clearly specifying what that GTF file will annotate. (It's kind of like saying "I need an XML file" or "I need an Excel spreadsheet" without indicating what the contents of those files would be.) The things you need to clarify are:

What will the intervals in my GTF file represent?
What precise sequences will serve as the reference and provide interval coordinates for the GTF file?
How are those intervals to be determined?

It sounds like 3. is a big part of your question. If I'm reading things correctly, you want to annotate a particular class of genes on an influenza genome. Is that correct? Did you perform de novo sequencing of influenza viruses, and now hope to characterize variation in their HA genes? Stating what you hope to get out of your analysis would be helpful in finding a good way forward.


how to sort and join two files based on First Column Id

You should check left_join() method in dplyr package. In this way you do not need to sort the data frames, as the function automatically take care of this.
For example: merged_df <- left_join(otu, taxa, by = "ID").


How to construct the phylogenetic tree of the "High Recombination" genus streptomyces?

You need bootstrap support using a model-based tree building algorhithm, via maximum likelihood (a few people use Bayes). The file format is relaxed phylip format (please submit a separate question if you have difficulties here - its a bit tricky).

I use standard RAxML here, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

raxmlHPC -m GTRGAMMA -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,

raxmlHPC -b 459 -m GTRGAMMA -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

This tree will bootstrap for 500 replications, however to start I would use 100 replications.

Make a consensus tree of the bootstraps,

raxmlHPC -f b -m GTRGAMMA -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)

The cluster you have access to is fine for the calculation, it will take aound 24 hours to complete one bootstrap calculation for one data set and obviously you need to parallelise you calculation across the 22 contigs.

Viewing the tree, FigTree (for Mac OSX) is easy.

Root

Rooting can be complicated because I don't really know your bacteria.

The recombination issue is more complicated, but I would construct 22 trees from your contigs and assess them for congruence. Panmixia is a concerned, which means too much recombination


Double-counting coverage of overlapped read pairs

samtools depth is a simplified version of samtools mpileup, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use samtools mpileup instead:

$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

I haven't shown the entirety of the output and you'd want to use the -aa flag to ensure that all bases are returned.

Note that samtools mpileup is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to -Q 1 (the default is -Q 13, which you'd want to change).

As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use bamPEFragmentSize from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.


Processed spliced and unspliced count matrices for existing scRNA-seq atlases

I am not sure if there are preprocessed datasets available, however, Kallisto + Bustools was designed specifically for the rapid quantification of spliced and unspliced transcripts.


Tumor sample and heterozygous SNVs: what is the genotype of the normal cells?

Your initial information/image and actual question seem pretty unrelated to me, but I'll do what I can.

If you have a heterozygous SNV in a malignant population of cells, it really depends on when the mutation occurred for you to know the genotype of wildtype cells. It could have been there since birth (a germline mutation passed down from the parents) or have been acquired throughout the person's life and clonally expanded (somatic mutation).

Unfortunately, you have no way to differentiate between these two types of SNVs just by sequencing the tumor cells. If you want to know true somatic mutations that are found only in the tumor cells, you need to sequence a matched normal control from the same patient (preferably from the same cell type). This is why most large cancer sequencing studies include not only tumor samples, but normal cells from the same patient - see TCGA or ICGC for examples.


NCBI retrieve genome assembly with splicing or "modifications"

The page you link to is a GenBank Flat file, the default format used by GenBank. So, to get that format, which includes the join lines, use:

./efetch -db nuccore -id CP003820.1 -format gb

The join() lines give the different ways in which a sequence can be built from this gene. For example, one mRNA sequence produced by the locus you mention (CP003820) is:

  mRNA            complement(join(100..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..>5645))

The corresponding coding sequence is for the mRNA above is:

 CDS             complement(join(168..1497,1556..1767,1823..2274,
                 2322..2782,2846..3126,3206..3890,3958..5263,5322..5422,
                 5494..5645))

These aren't necessarily splicing "modifications", and aren't "variants", they are simply the way that the exons are joined to make an mRNA.

Note that since CP003820.1 is "Cryptococcus neoformans var. grubii H99 chromosome 1, complete sequence", there are multiple genes described in there and each will have its own (or multiple, if there is alternative splicing) join statements. Specifically, I see 1014 genes described in the file:

$ grep -Pc '^\s+gene\s+' CP003820.gb
1014

950 mRNAs:

$ grep -c 'mRNA.*join(' CP003820.gb
950

And 955 proteins:

$ grep -Pc 'protein_id=' CP003820.gb
955


How do I change the identity of a sample from spatial 'Visium' data preprocessing in Seurat v3?

You can re-assign orig.ident after creating the object. For example:
brain1 <-Load10X_Spatial(data.dir = "~/xxx", filename = "filtered_feature_bc_matrix.h5", assay = "Spatial", slice = "slice1", filter.matrix = TRUE, to.upper = FALSE)
brain1$orig.ident <- "Visium"


Retrieve overlap widths with Bioconductor

Whilst intersect is sort of doing what you want, it merges adjacent overlaps together which is not what you want I suppose. There doesn't seem to be dedicated Grange function to do this out of the box, but you could go through all overlapping pairs resulting from findOverlap and calculate it's widths using pintersect:

punion, pintersect, psetdiff and pgap are generic functions that compute the element-wise (aka "parallel") union, intersection, (asymmetric!) difference and gap between each element in x and its corresponding element in y.

"pintersect"(x, y, resolve.empty=c("none", "max.start", "start.x"), ...) # x, y   IRanges objects.

hits <- findOverlaps(gr1, gr2)

width(pintersect(gr1[queryHits(hits)],
                 gr2[subjectHits(hits)]))
[1]  2  6  4  6 10


How can we find the distance between all residues in a PDB file?

If you need to process multiple files, you could use Biopython to parse a PDB structure.

from Bio.PDB import PDBParser

# create parser
parser = PDBParser()

# read structure from file
structure = parser.get_structure('PHA-L', '1fat.pdb')

model = structure[0]
chain = model['A']

# this example uses only the first residue of a single chain.
# it is easy to extend this to multiple chains and residues.
for residue1 in chain:
    for residue2 in chain:
        if residue1 != residue2:
            # compute distance between CA atoms
            try:
                distance = residue1['CA'] - residue2['CA']
            except KeyError:
                ## no CA atom, e.g. for H_NAG
                continue
            if distance < 6:
                print(residue1, residue2, distance)
        # stop after first residue
        break

If you need to look at one structure, using a viewer perhaps would be easier. You could try PyMOL: (how to measure distance). There are other PDB viewers, some of which can work even through a browser.


A new paper suggests the Corona Virus has "Uncanny similarity of unique inserts in the 2019-nCoV spike protein to HIV-1" - What does this mean?

UPDATE: The article has now been withdrawn with the following note:

  This paper has been withdrawn by its authors. They intend to revise it
  in response to comments received from the research community on their
  technical approach and their interpretation of the results. If you
  have any questions, please contact the corresponding author.

This is very odd, and will require a rigorous investigation, but my initial reaction is one of scepticism. Considering just the 1st insert, the insert sequence is GTNGTKR, short at just 7 amino acids. A simple BLASTP vs NR did not find perfect matches to HIV sequences, but did reveal 100% identity across the full sequence length to >50 other short protein sequences, which could be spurious 'chance' hits of course. Many of the top 100 hits are against eukaryotic sequences, for example one is against Pristionchus pacificus, a type of nematode worm (see alignment below). Someone needs to do a proper peer review of this preprint before any conclusions are drawn from it.

>tank-1 [Pristionchus pacificus]
Sequence ID: PDM74036.1 Length: 2481
Range 1: 1474 to 1480

Score:24.0 bits(49), Expect:2477,
Method:,
Identities:7/7(100%), Positives:7/7(100%), Gaps:0/7(0%)

Query  1     GTNGTKR  7
             GTNGTKR
Sbjct  1474  GTNGTKR  1480

UPDATE 1: Things are moving fast. There are now 10 comments under the preprint, all of them critical of the idea that these inserts are meaningfully similar to HIV-1 sequences for the reason I outlined above (i.e. short sequences with many hits against organisms from across the tree of life) and additional points, such as small insertions being quite normal evolution for RNA viruses. There are also similar critiques on twitter and biorxiv has added the following disclaimer header to their website:

  bioRxiv is receiving many new papers on coronavirus 2019-nCoV.   A
  reminder: these are preliminary reports that have not been
  peer-reviewed. They should not be regarded as conclusive, guide
  clinical practice/health-related behavior, or be reported in news
  media as established information.

I have never seen such rapid (and unanimous) post-publication peer review before. It is therefore already (just 1 day after publication) fairly clear that at least some claims in this study are completely false and there is no good evidence to support the claim that the 2019-nCoV have acquired sequences from HIV or evidence to suggest the virus is engineered.


Extracting modules from correlation and then module clustering

Co-expression in R is doable via gglot2 at the very least... I see if I can find the code.

Essentially you are performing a co-variance analysis, however Pearson' at 0.4 is a LOW cut-off and it should be 0.7.

I've got a lot on right now, so it will take a week or so before I foward the relevant code. R will have other methods outside ggplot2 however.

I was thinking ggcorr() and ggpairs() which both in the GGally package of R-cran. The output is very cute, well its ggplot2 which is trendy.

library (ggplot2)
library (GGally)
library (dplyr)

df <- select(mydata, 1:5)
cor(df) # raw output
ggcorr(df)

output: heatmap similar to Fig. a

ggpairs(df)

output: very fancy matrix plot Its pairwise scatterplots, diagrams and numbers, which is not represented in Fig. a to c, but its a very fancy version of Fig. d ... Its simplistic analysis, but its a very cute plot.

ggcorr uses pairwise (obviously) Pearson by default, which is the most versatile.

I can definitely say that 0.7 for Pearson for a straight correlation (not a covariance analysis) is conservative low and as a reviewer I would immediately complain about anything lower. A heatmap is different depending on its interpretation.

The multivariate analysis of Pearson pairwise correlation, which is what you are requesting is a different question. UPGMA style analysis directly on the heatmap, is fine if there is "one heatspot", but otherwise has limitations IMO.


Unfortunately, PDB numbering can be one of different schemes. The main three:

Uniprot numbering (canonical isoform)
Construct numbering with or without tag removal
Start from 1 for what has density

However, the are many many corner cases. In Proteopedia there is a lovely article about weird PDB numbering.
The best way to find the details about what part of the sequence is covered is to use the Swissmodel API. It's actually better than either the PDBe and RCSB PDB APIs. Swissmodel is a homology modelling service that maintains close threaded models (>30% ID), but also has PDB data. Here is a Python3 method I use, (provider='pdb') is what is required here).
def get_data(uniprot: str, provider='swissmodel') -> List[dict]:
        """
        Gets data from Swissmodel. Convert each entry with convert_SM_entry
        """
    assert provider in ('swissmodel', 'pdb'), 'Provider has to be swissmodel or pdb'
        url = f'https://swissmodel.expasy.org/repository/uniprot/{uniprot.strip()}.json?provider={provider.strip()}'
        data = requests.get(url).json()
        return data['result']['structures']  # [0] # from to coordinates description

The returned dictionary contains a start and end of the uniprot <--> PDB. It does not however, say where the missing density residues are. But in deposited structures the PDB numbering skips where these residues ought to be... So might not be an important problem.
Another database of interest is SIFTS.
To make a multiple sequence alignment of a non-deposited PDB, you need to extract the sequence and do the alignment here is a snippet using pymol module (from conda):
import pymol2

pdbfile = '' #your pdb file name..
selector = 'chain A' # or whatever chain you want
with pymol2.PyMOL() as pymol:
    pymol.cmd.load(pdbfile)
    fasta = pymol.cmd.get_fastastr(selector)

seq = ''.join(fasta.replace('?','').split('\n')[1:]) # remove header and non AA residues

Biopython has a dozen sequence alignment tools. But the original one still works fine
from Bio import pairwise2
from Bio.Seq import Seq

alignment = pairwise2.align.globalxx(''.join(seq).replace('?',''), "ACQQG")[0]

alignments is a tuple of the MSA sequences.
If you want to do the latter without coding. Extract the sequence from the PDB using PyMOL the normal application and load the pdb and then print cmd.get_fastastr('chain A') in the command line part of the GUI. Copy this. Go to Muscle or other onliner aligner and align!


How to demultiplex a mix of single-indexed and dual-indexed samples

I'm late to the party but simply run cellranger mkfastq twice, with different arguments for the mask and with --filter-dual-index for the double-indexed samples. In a second step you then have to disentangle the dual A+B reads from single-index A reads based on the readID's that are found in both of the demultiplexings. This is cumbersome and time+diskspace consuming, but it can be done. I would use some Unix tools for that but there are probably smarter ways.
EDIT: I mistakenly said bcl2fastq, this must be cellranger mkfastq. (thanks AmadeusDrZaius)


Is there a way to assemble contigs starting from a specific sequence?

This is one of the primary use cases for which spacegraphcats (preprint and code) was designed. The "neighborhood queries" discussed in the paper sound particularly relevant. I don't have any personal experience with spacegraphcats, but the run guide provides some examples of how to index the complete data set and how to query sequences of interest.


1000 genomes missing SNPs on chr X in GRCh38

As I told you on 1000 genomes helpdesk:
Because of the ploidy settings we used for our variant calling, we only called variants in the the pseudoautosomal region on the X chromosome. We are working on settings to correctly analyse X variants in both males and females and plan to analyse X properly soon.


Can't locate object method "new" via package Bio::Tools::Primer3

I managed to figure out the answer to my question. I installed Bio::Tools::Primer3 by typing the following into terminal: 'sudo cpanm Bio::Tools::Primer3 --force' to install it. It only works when I used 'sudo' and 'force'


Tag density plot form chip or atac seq data

You can use Homer to retrieve tag densities, quoting from their manual:

  Calculating ChIP-Seq Tag Densities across different experiments

  annotatePeaks.pl is [a] useful program for cross-referencing data from
  multiple experiments.  In order to count the number of tags from
  different sequencing experiments, you must first create tag
  directories for each of these experiments.  Once created, tag counts
  from these directories in the vicinity of your peaks can be added by
  specifying "-d   ...".  You can
  specify as many tag directories as you like.  Tag totals for each
  directory will be placed in new columns starting on column 18.  For
  example:

  annotatePeaks.pl pu1peaks.txt mm8 -size 400 -d Macrophage-PU.1/
  Bcell-PU.1/ > output.txt


Functional Annotation vs Functional Enrichment--which one should I use for Network Analysis?

The answer to your main doubt about "annotation vs. enrichment" is that you need both, but you already have one annotation.
You can only perform a functional enrichment analysis after you have a functional annotation. In your question, you state that you have KEGG annotation for your targets. That is it, knowing which target "belongs" to each KEGG module is an example of functional annotation of your targets. Several tools can annotate genes, proteins, etc. Knowing which protein is related to each drug or disease is also a functional annotation.
Now, knowing which targets are significantly important can be achieved by functional enrichment analysis.
For the network, both information can be used together. For example, you could link all elements to their functional annotation and use the available enrichment score as a quantitative parameter of importance for calculating correlations.


Get list of urls of GSM data set of a GSE set

You can use Entrez Direct for this as follows:
esearch -db gds -query 'GSE104279' \
  | esummary \
  | xtract -pattern DocumentSummary \
    -if 'entryType' -equals 'GSM' \
    -def 'NA' -element Accession title summary FTPLink

This will return a table with data similar to this:
GSM2580330      18_Z13_2_d5_Zika        cortical organoids_Zika_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580330/
GSM2580329      17_Z13_2_d5_Control     cortical organoids_mock_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580329/

To download sequence reads, you should follow links to SRA. Using Entrez Direct you can do this as follows:
esearch -db gds -query 'GSE104279' \
  | elink -target sra \
  | efetch -format runinfo

This will return a comma-delimited table containing SRA identifiers and an FTP path to the SRA data. These FTP paths won't lead you to FASTQ files though. You can pass the SRA run identifiers of the format SRR### (or ERR### or DRR###) to fastq-dump or fasterq-dump tools from the SRAToolkit to download data in FASTQ format.


How can I obtain FTP links to studies in ENA?

Here are some command lines I used for that purpose in bash.
Simply prepare a text file containing each accession number (SRR/ERR) you want and create a for loop. Here I used prozilla to speed up downloads but you may use wget either.
for index in $(cat list_of_accessions) ; do

if [ ${#index} -eq 9 ]; then
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_1.fastq.gz
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_2.fastq.gz
            else
            if [ ${#index} -eq 10 ]; then
                    proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_1.fastq.gz
                            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_2.fastq.gz
                    else
                    if [ ${#index} -eq 11 ]; then
                     proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_1.fastq.gz
                                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_2.fastq.gz
                    else
                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_1.fastq.gz
                                                proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_2.fastq.gz
                        fi
                fi

            fi
done
```


Are there any tools for discovering tandem repeats around a variant in a VCF file?

After researching a bit, I found a couple of tools that take as input BAM files. I did not test them yet:

GangSTR
lobSTR
hipSTR
Tandem Repeat Finder (thanks to Michael G.'s answer)

If I find more tools, I will update this post. After testing them, I will come back to share my feedback.


Restriction enzymes with multiple recognition sequences

I would honestly just remove the non-homologous cutting site enzymes from the database. Firstly, these are biology students and require a simple picture to grasp the concept of coding via regex pattern matching. Secondly, in a wet-lab nobody uses these enzymes, I was not even aware the "8" existed until your post. In a teaching lab the best motto is "keep it simple".


How to input data for DESeq2 from individual HTSeq count?

In case still needed, here is a piece of code that I have used for the same purpose (includes a lot of sanity checks). This is highly likely an overkill as each count file would probably have the same rows (as a result of being generated from the same workflow) and thus can easily be merged column-wise.

count_file_dir <- "../input/count_files/"
count_file_names <- list.files(count_file_dir)
count_files <- paste0(count_file_dir, count_file_names)

# read and merge count files, while removing last 5 lines that correspond to summary stats
count_table_list <- mclapply(count_files, function(x) {fread(x)[1:(.N-5)]})
names(count_table_list) <- sub(".count", "", count_file_names)

# check the dimensions of count files
t(sapply(count_table_list, dim))

# check if the count files have the same gene ordering
# compare gene names of the first file with all
gene_names <- sapply(count_table_list, function(x) x$V1)
as.data.frame(apply(gene_names, 2, function(x) all(x == gene_names[,1])))

# add file name as a unique column name for gene counts
for(i in names(count_table_list)) {
  names(count_table_list[[i]])[2] <- i
}

# merge read counts and convert to a matrix
count_matrix <- Reduce(merge, count_table_list)
count_matrix_rownames <- count_matrix$V1
count_matrix$V1 <- NULL
count_matrix <- as.matrix(count_matrix)
rownames(count_matrix) <- count_matrix_rownames
colnames(count_matrix) <- sub(".count", "", colnames(count_matrix))

dim(count_matrix)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                              colData = ...,
                              design = ...)


Why isn't there a standard method to convert GFA to JSON?

Short answer: no, there isn't a "standard" method, and there won't be.

Long answer: like TAB-delimited and XML, JSON alone is not a specific format. You have to define a schema to give meanings to data. However, unlike TAB-delimited "format" and XML, there is not an official way to define JSON schema. There have been various attempts but all of them are complicated and none of them are widely adopted. As you can't write down a JSON schema for GFA in the first place, there won't be a "standard" way to convert GFA to JSON. The best you can do is to roll out your own in-house solution, unfortunately.


Seurat define manually a cluster and find markers

It would have been better if you have shared your code as well as the error messages you got format these as such.

As far as I can understand, your problem is that you failed to make your cluster as default ids. For this you can use:

seurat_obj <- SetIdent(seurat_obj, id=seurat_obj@meta.data$column_of_interest).

The cluster identities are stored in the @meta.data slot of your object.


Storing FindAllMarkers results in Seurat object

You can stash anything you like in the misc slot (present in both v2 and v3 Seurat objects). I often use it for storing marker information to help with organization.


Retrieve abstract/summary for NCBI bookshelf via Entrez

Entrez requires the unique identifiers(UID) for fetching related info. The id you are using in the query is RID. I guess that why you are getting <ERROR>UID=1440: cannot get document summary</ERROR>.

List of example uids using the term "ncbi+blast":

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=books&term=ncbi+blast

Summary of the first item:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=books&id=4459608

Using the Entrez python package:

from Bio import Entrez
Entrez.email = "youremail@domain.com"
handle = Entrez.esummary(db="books", id="4459608,4456634", retmode="xml")
records = Entrez.parse(handle)
for record in records:
     print(record['Title'])
handle.close()

I'm not sure about the availability of book abstracts.


differentially expressed intron analysis?

Typical tools to use for this would be the DEXseq packages (as done here) or rMATS, which will calculate intron retention rates.

I would generally suggest not summing intronic counts by gene and using that in DESeq2 (I assume that's what you meant by DESeq, since there's no reason to use the original DESeq for anything). What you'd be getting if you did this is less expressed introns and more sample-specific bias in enriching for mature versus pre-mRNA.


Find intervals that genes fall within their range

The GenomicRanges package is a great place to start, something like this:

library(GenomicRanges)

# make a GenomicRanges set corresponding to the genes
genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom")

# turn your windows into a data frame and make a GenomicRanges set from that
windows = ...
windows.gr = makeGRangesFromDataFrame(df = windows ...

# find overlaps between the two sets
olaps = findOverlaps(genes.gr, windows.gr)

# output relevant fields
cat(paste(genes[queryHits(olaps),"gene_id"], windows[subjectHits(olaps),"window_id"]))


WGCNA error codes in Network Construction & moduel detection

Looks like a bug in the code which I will try to hunt down and fix. In the meantime, I would suggest playing with maxBlockSize argument to blockwiseModules. Try increasing it as much as your available RAM allows (see the paragraph "A second word of caution concerning block size" on page 6 of WGCNA tutorial I, section 2c, for some guidance on how to set maxBlockSize). Hope this helps.