Microarray Meta-Analysis and Cross-Platform

In short the answer is no because you have to standardise between the arrays. Thus both arrays would need to run exactly the same control sample and you would then perform statistical standardisation (its a formal term), which is a varience measure where you look at the variation around the 'mean' for all samples in the array and here the 'mean' is the standard sample. If that can be done for both arrays at that point they can be combined.
The problem is intensity per run varies and you have to control against this

In this scenario you  can certainly combine samples by assessing whether the control samples between runs approximates to identity. What you could use is a uniform distribution to assess this, variance must MUCH smaller than the mean for the STANDARD sample(s). I don't really do uniform distribution statistics, so I'd need to research this, because 'uniform' means identity to standard deviations against the uniform distribution is a bit meaningless. However again the uniform distribution is denoted by a very variance in comparison to the mean. The statistic of dividing mean by standard deviation/variance from memory is called the Index of Association (I could be wrong). If this value is at an extreme you're good to go <<<1, if it is remotely close to 1 you're in trouble, if it is >1 you need to use standardisation.
Presumably you have several standards to ensure a mean is meaningful. A T-test and non-parametric equivalents is not stringent enough. If there is a difference assess a correction factor and rerun the uniform distribution check. The same correction will be applied to all samples.
The undisputed method to combine is standardisation, but accept in this scenario the differences between runs is potentially very small, so here what you are attempting is essentially to combine raw data or at least keep one of the runs as pure raw data.


run a tool with python>=3.6 on galaxy-project instance that runs on python2.7

Each Galaxy tool uses its own (typically conda-based) environment, which will normally not use the same python as Galaxy itself (note that Galaxy can run under python 3). An example of this is deepTools, which is python 3 only but will happily run under a python 2 Galaxy since conda is being used to create a python 3 environment for it. It's best to distribute your tool via bioconda and then Galaxy will be easily able to install it.


Accuracy of the original human DNA datasets sequenced by Human Genome Project?

The HGP developed the first "reference" human genome - a genome that other genomes could be compared to, and was actually a composite of multiple human genome sequences.

The standard human reference genome is actually continually updated with major and minor revisions, a bit like software. The latest major version is called GRCh38, was released in 2013, and has since had a number of minor updates.

  Are the datasets provided by HGP still accurate?

Yes, in a sense, but we certainly have better information now. One way to measure the quality of the assembly is that the initial release from the HGP had hundreds of thousands of gaps - sequences that could not be resolved (this often occurs because of repetitive sequences). The newest reference genome has less than 500 gaps.


Is there anything similar to GSEA for locus-based (instead of of gene-based) data?

I tinkered with a program GSEA-SNP quite a few years back, which claims that it does a similar ranking procedure with SNPs. It carries out its procedure by first linking SNPs with genes, then running an algorithm similar to GSEA. There's a bit more detail in the paper.

Unfortunately this is a space I've got a bit further away from in recent years, so I don't know if anyone has carried out the legwork to generalise the algorithm for sets of any type where each item within each set fits a continuous (but not necessarily normal) distribution.


How can I classify the 3 clades(S, G, V) of the coronavirus without using protein data?

It seems to be explained right there in the image you posted:

So, the three strains were classified based on three specific variants:

Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the  gene's protein product   to a valine (V).


biological fluorescence supervised machine learing

Whether to quantify things at all with depend on the technique you want to use. An end-to-end CNN (convolutional neural network) for example would just be fed the 3D image for training and predictions. However the major downside to CNNs is that you need a LOT of data to train them, though perhaps you can get lucky and find an existing CNN that you can use and just retrain the last couple layers.

Other than that, you're going to need to come up with as many (hopefully useful) metrics as you can. There are then a lot of techniques one could try, ranging from simple things like regression to random forests to SVMs, among many others. In general I think simple regression and random forests are the easiest to get started with.


Is it possible to use spike-ins with droplet technologies?

I believe it is mostly a cost issue as spike-ins are "expensive" to use given the large number of cells that are/can be sequenced via droplet based methods (as compared to plate-based ones). The 10x link that you provided actually states that this is possible with 10x and Klein et al. used ERCC spike-ins in 953 droplets (Cell. 2015 May 21;161(5):1187-1201. doi: 10.1016/j.cell.2015.04.044.).


I was able to reproduce your problem, even in a blank environment.

$ conda create -n test
$ conda activate test
$ conda install bioconda::blast=2.7.1
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - bioconda::blast=2.7.1
  - boost=1.64
  - bioconda::blast=2.7.1
  - gnutls
  - bioconda::blast=2.7.1
  - nettle[version='3.3|3.3.*']
  - bioconda::blast=2.7.1

Note that boost, gnutls, and nettle are listed as requirements here.

If you do a search for those packages, you'll see that boost is available in a default channel, whereas gnutls and nettle are only in the conda-forge channel, which is not a default.
Adding those packages to the install command resolved this issue, for me.

$ conda install bioconda::blast=2.7.1 conda-forge::gnutls conda-forge::nettle --dry-run
Solving environment: done

## Package Plan ##

environment location: /Users/hawleyj/anaconda3/envs/test

added / updated specs:
    - bioconda::blast=2.7.1
    - conda-forge::gnutls
    - conda-forge::nettle
...

As mentioned in the comment below, to avoid this issue entirely, before using Bioconda packages make sure you run the proper channel configuration first.

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge


A tool or webserver for building PSSM matrix

If you take a look at my answer in this BioStars post, you can generate a PSSM using AlignIO in Biopython:


Why is the UniProt REST API returning multiple results, when I am only providing one ID to be checked?

it is because you do a text search and the "word" Q9UJL9 appears in multiple entries. If you want an id search, specifically say so and for that the query syntax is id:Q9UJL9. You will want to follow a redirect, and there are some rare cases where that might still lead to more than one entry if the id was made secondary and attached to multiple entries. e.g. P00001.


Histogram of different classes of variables in one column

This isn't a bioinformatics question, but it's quicker to answer than to close.

df1 <- read.table(text = "
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7", na.strings = "")

library(ggplot2)

ggplot(df1, aes(x = V1, y = V2)) +
  geom_bar(stat = "identity")


How do I force Cromwell / WDL to put input files in the same directory?

My hacktastic fix for this is to use basename to get the get the name of each file, then do
refName=`basename "${Ref}"`
ln -s "${Ref}" "$refName"

for each file.  This puts aliases to each file in the executions dir, and then I pass "refName" to GATK rather than "${Ref}".
Not a good solution, but a solution


Your shell script is fine, and works as expected. The reason it hangs is that you aren't telling PyMol to exit after running your script, so it's still running and waiting for the next command. Just add a quit at the end and it should work:

#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
select prot, 1OYT and not ligand,
save lig.pdb, ligand
save prot.pdb, prot
quit
' > scriptt.pml
pymol scriptt.pml
obabel prot.pdb -xr -O prot.pdbqt
obabel lig.pdb -O lig.pdbqt
echo 'babel conversion is finished'


There is a recent paper that attempts to do this:

  ISOWN: accurate somatic mutation identification in the absence of
  normal tissue controls.

  In this work, we describe the development, implementation, and
  validation of ISOWN, an accurate algorithm for predicting somatic
  mutations in cancer tissues in the absence of matching normal tissues.
  ISOWN is available as Open Source under Apache License 2.0 from
  https://github.com/ikalatskaya/ISOWN


How to load output from the kegg API in biopython into a pandas dataframe?

solution from @sören:
I came up with this solution for my problem:
from Bio.KEGG import REST as kegg

def _get_kegg(kegg_id):
    kegg_output = kegg.kegg_get(kegg_id).read()
    results = {}
    for line in kegg_output.split('\n'):
        splits = line.split()
        if not line.startswith(' '):
            if len(splits) > 0:
                key = splits[0]
                value = ' '.join(splits[1:])
                results[key] = value
        else:
            results[key] += ' '.join(splits)
    return pd.DataFrame(results, index=[kegg_id])

_get_kegg_v = np.vectorize(_get_kegg)

def get_kegg_info(kegg_ids):
    if isinstance(kegg_ids, str):
        kegg_ids = [kegg_ids]
    return pd.concat(_get_kegg_v(kegg_ids), sort=False)


Clarification about how mapping to a reference genome works in RNA seq

You can't always know exactly where every read originated, even if you use long reads to sequence full-length transcripts (odds are there's at least one transcript that's duplicated in the genome). This is also not exactly a problem, such reads can either be discarded or reported as multimapping.


Can ortholog genes have different names among different species?

Gene naming is pretty subjective and even a single gene can have many names. If you're lucky orthologs will be named the same in different lineages, but you can't rely on it. Better to use gene ids (e.g. from Ensembl) and an ortholog db to map the orthologs. However, if you're looking at different strains, this is more tricky since most db don't go below the species level


How to split multifasta based on partial fasta header

The SeqBuddy --pull_records command can grab sequences based on a regular expression against record IDs and/or metadata. Given your example:

$: seqbuddy.py sequences.fa --pull_records "Clostridium botulinum.*BKT028387" full > Clostridium_botulinum_BKT028387.fas
$: seqbuddy.py sequences.fa --pull_records "Peptoanaerobacter stomatis.*OBRC8" full > Peptoanaerobacter_stomatis_OBRC8.fas


What are the right parameters to trim a small RNA transcriptome with trimmomatic?

The command where you trim with adapters and by quality is perfectly fine. That FastQC isn't perfectly happy is expected. It's a tool made for whole genome sequencing QC, you SHOULD see a number of "fail"s with any RNA-seq protocol. Steps that should always fail in RNA-seq:

per-base sequence content (there's "random priming" performed that isn't completely random, so the first ~8 bases show a typical bias pattern).
per-sequence GC content (particularly in sRNA-seq you're enriching for a subset of genes, they won't have some nice GC distribution and that's perfectly OK).
Sequence duplication level (if this doesn't fail then everything is uniformly expressed...which would be quite unlikely)
Overrepresented sequences (as above, though this will be worse for sRNA-seq since reads will often cover the entire transcript, so there's even less sequence diversity)


How to run a one-time 'setup' step for amrfinder in snakemake using conda

There are two approaches you can take:

There's a --conda-create-envs-only option that you can use to exit after creating the envs. You should get information about where the envs are so you can run amrfinder -u.
Add a new rule called "amr_setup" or something like that that runs:

.

amrfinder -u
touch amrfinder.setup

Then put amrfinder.setup into your snakemake workflow as an input file into the amr_finder_plus rule. Either of these will work


How to translate amino acid sequences to Nucleotide sequences

Although there is not a unique nucleotide sequence that translates to a given protein, one can list all the possible DNA sequences that do translate to that protein.

An online tool that does just that is Backtranambig, from EMBOSS. It produces a DNA sequence representing all the nucleotide sequences matching the input protein, using IUPAC ambiguity codes.


Why do NEBNext indexing primers have sequence between the p5 oligo and index?

Question 1: The additional sequences is needed because that is complimentary to the P5 sequence anchored to the flow cell. It is also the site for priming the Index 2 read on a MiSeq. The additional sequence is also needed for the read 1 primer in the cartridge to anneal to the correct place for the molecules.

Question 2: The "GT" sequence you're referring to is necessary for the Index 1 and Read 2 primer to anneal to the correct place.


Are mitochondrial genes to exclude in scRNA-seq such as ribosomal genes?

According to Ilicic et al. (2016), on upregulation of mtRNA in broken cells:

  There is an extensive literature on the relationship between mtDNA, mitochondrially localized proteins, and cell death [34, 35]. However, upregulation of RNA levels of mtDNA in broken cells suggests losses in cytoplasmic content. In a situation where cell membrane is broken, cytoplasmic RNA will be lost, but RNAs enclosed in the mitochondria will be retained [...]

According to this, I think that mtRNA levels relative to endogenous RNA can be used as a control for low quality (broken) cells, in a similar way to ERCC spike-ins.

Ilicic T, Kim JK, Kolodziejczyk AA, Bagger FO, McCarthy DJ, Marioni JC, Teichmann SA. Classification of low quality cells from single-cell RNA-seq data. Genome biology. 2016 Feb 17;17(1):29.


You don't have blastp in your path. Blast is either not installed on your system or it's just not been added to path right.


EDIT after OP's update:

Try this Perl one-liner:

perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file

Here, the command line flags are:
-e: tells the perl interpreter to look for the code inline, rather than in a file with the script,
-p: loop over the input one line at a time, execute the code in the one-liner for each iteration of the loop, and print the result into STDOUT.
The code has regex replacement (see perlre for more details):
s{>([^_]+)_.*_}{>RNA${i}#${1}/}:
replace an underscore, followed by zero or more occurrences (* quantifier) of any character (.), followed by an underscore, with slash (/).
Keep > and the original header part before the underscore, captured by parens into $1 variable.
Add string RNA${i}, where $i is the counter of the fasta headers. The counter is incremented by 1 when the substitution occurs, that is, when the next header occurs.
The variables are used like so in the substitution: ${i}, ${1} to prevent wrong interpolation.
I am also using s{}{} as delimiters for readability (as in most cases) instead of the default s///.

EXAMPLE:

echo ">LTR_retrotransposon100_Gypsy?\nACGT\n>LTR_retrotransposon10011_ERVK\nTCGA" | perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}'

Prints:

>RNA1#LTR/Gypsy?
ACGT
>RNA2#LTR/ERVK
TCGA


Is there an established method for comparing expression of groups of genes (gene sets)?

You could apply the GSVA transformation to later compare them. You would need to give a step back and start from the expression values not log2FC, but I think it would be better than using the log2FC directly.
If you need to use them you could compare the gene set enrichment value of these two sets on the same comparison. I would use fgsea with the log2FC and the two sets of genes. If the normalized enrichment value (NES) is different then you could explain that these sets expression are indeed different.


Writing a perl script to holding information for two genes

Perhaps this small test script will help demonstrate some of the principles:

#!/usr/bin/env perl

use strict;
use warnings;
use Data::Dumper;

my @arr;

# define the current line
my $line = "foo\tbar\tbaz\n";

# remove the trailing newline character
chomp $line;

# put the line's elements into an anonymous array
my ($a, $b, $c) = split('\t', $line);

# insert key-value pairs into anonymous hash
my $hash = {a => $a, b => $b, c => $c};

# add the hash table to the array
push @arr, $hash;

# print the array
print Dumper \@arr;

Run this script and see what it prints out.


I am unaware of any "official" or gold-standard UTR annotations in S. cerevisiae.

One option is to use the annotations from the TIF-Seq publication (Pelechano et al. 2013).

The GSE39128_tsedall.txt.gz file contains the major isoforms identified. It would be up to you to computational associate each transcript with a given gene. It is also up to you to determine which major isoforms are present in your dataset. Furthermore based on my own experience there is often a significant number of discrepancies between what I see in an NGS dataset (e.g. RNA-Seq) and what is called as a possible isoform in these annotations. However, a lot of this could be related to the inherent technological limitations of the method.

If you look at the SGD genome browser they provide some additional UTR annotations. Derived from Nagalakshmi et al. 2008.

You can find the relevant annotations for those UTR's in their supplemental data section. The file is 1158441_tables_s2_to_s6.zip and inside you will find TableS4 which seems to contain what you seek.

There might be more resources but these are the two main ones I'm aware of.


How to force pilon to correct a reference with low coverage mapped reads?

Pilon will not consider fragments that are not proper pairs if the --frags option is used to pass the alignment .bam file. The reads mapped in the above IGV screenshot are unpaired.

Instead, if you want Pilon to consider all reads that map as legitimate and usable for correction, pass the .bam file with the --unpaired option instead.

To above corrected example would be:

java -Xmx100G -jar {params.pilon}\
  --genome {input.ref}\
  --unpaired {input.bam}\
  --output {params.outname}\
  --changes --fix all,amb,breaks,bases\
  --threads {threads}\
  --mindepth 1


How to plot circular chromosomes with Gviz?

The primary author of Gviz confirmed to me via email that this is not possible:

  Gviz will only give you the traditional linear tracks. There are a
  bunch of solutions available for visualizing circular genomes, most of
  which are more interactive in nature and not implemented in the rather
  static R graphics context. GGplot2 has some support for circular axes,
  too. If you want to stay within the R environment I’d start there.
  Ggbio is another option providing a somewhat higher-level interface.


How to fix 5'UTR annotation with RNA-seq data?

Genes, particularly transcribed genes in Eukaryotes, can have alternate start and stop sites (and splicing).  Alternate transcriptional start sites may be found in different tissues or at different developmental periods due to alternate transcription factor binding sites being used.  The only reliable way to determine the transcriptional start and stop sites is by read coverage.
Genome annotations are generally the product of bioinformatic pipelines and are mostly predictions based upon computational models.  Different gene-finders may call different boundaries; for example, in eukaryotic genomes, an exon may be assigned to the wrong up/down-stream gene.  Arabadopsis, as model organism, has had it's gene boundaries considerably refined/corrected/validated, but even so, there could still exist errors in the annotation or, more likely, the boundaries indicated are the most common start/stop sites (see point 1 above).

In short, your UTRs are whatever your experimental data (reads) show it to be.  Don't worry too much about it not matching up to the standard annotation boundaries.


Parsing the json file from gene ontology result into dataframe

I don't know how to do that in R but I think this is what you want:
I used http://json2table.com for view the data structure of the JSON file, then I used json_normalize to normalize the JSON data into a flat table using the key [result] to a list of DataFrames for each index in the JSON file. Then, I concatenated the list of DataFrames in a single one with df.concat()
import json
import pandas as pd

with open('analysis.json') as f:
    d = json.load(f)

dfs = []
for x, y in enumerate(d):
    data = pd.json_normalize(d[x]['result'])
    dfs.append(data)

df = pd.concat(dfs, ignore_index=True)

df.to_csv('analysis.csv', index=False)


bcl2fastq2 v2.20.0 converstion software on Ubuntu 18.04

Install it with conda.
conda install -c dranew bcl2fastq


Use a specific python.exe instance with PyMOL/Windows?

So it is really hard to make a given python instance run a library from a different instance. It is not just one dist-packages library folder, but multiple that you would need to link, and all that assuming they use the same version of python.
So the best bet is to install a new instance of pymol the python module.

The module can be installed in a python that uses conda with

conda install -c schrodinger pymol

This is the incentive version of PyMOL 2 and does not require a licence.
You installed rdkit, so I am assuming you used conda as conda install -c rdkit rdkit. But just in case you went through the pain of compiling rdkit for a given python instance without conda and want the same for PyMOL. There are guides of instructions required to compile open source. Such as https://pymolwiki.org/index.php/Windows_Install. On a barebone OS like a Pi's Raspsian, I had to do a lot more than claimed say sudo apt-get install python3-dev libglm-dev freeglut3-dev libglew-dev libpng12-dev libfreetype6-dev build-essential libxml++3.6-dev and link a few things manually. So just like for rdkit, conda make life a breeze.

PS. RDKit and PyMOL? Add PyRosetta and you just my favourite modules!


Plot to show the expression of genes between tumor and normal

You'll want to use long-form data for everything:

library(dplyr)
library(ggplot2)
logCPM$gene = row.names(logCPM)
d = logCPM %>% gather(Sample, logCPM, -gene)
d$group = c(rep("Tumor", 2), rep("Normal", 2))
geneOfInterest = d %>% filter(gene == 'RP11-351J23.1')
ggplot(geneOfInterest, aes(x=group, y=logCPM)) + geom_violin()

That's a simple example, you can tweak it to meet your needs.


Bootstrap analysis - why it is called bootstrap?

The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy of an estimate from random subsets of a population. Normally estimating something like the variance of the mean requires multiple independent samples, but bootstrapping allows you to perform estimates from a single population. This is essentially allowing the estimate to "pull itself up by its bootstraps", which is the only reasonable source for the term. See also the more general wikipedia article on bootstrapping outside of statistics.


Using Beagle 4.1 for phasing and IBD - chromosomes

OK, I ended up contacting one of the authors of the program, who kindly answered to my question:

  Either approach should work,

[i.e. running the analysis by chromosome vs joining the data together]

  but I suggest dividing by chromosome because that approach has had
  more testing.

I hope this might help other users, I saw that several questions regarding Beagle went unanswered on biostars and seqanswers. I guess it is a pretty specialized topic.

Best,

Max


How to convert species names into common names?

use the NCBI taxon dump under ftp://ftp.ncbi.nih.gov/pub/taxonomy

in taxdmp.zip you'll find all the names for a given NCBI taxon

$ grep -w ^9606 names.dmp
9606    |   Homo sapiens    |       |   scientific name |
9606    |   Homo sapiens Linnaeus, 1758 |       |   authority   |
9606    |   human   |       |   genbank common name |
9606    |   man |       |   common name |


I have found reason for not getting any result. I was using default adapter given by mirdeep* which gives me zero alignment with my data. This error pointed out by author of mirdeep* ( thanks to Jiyuan and colleen). After giving right adapter, I got the right result.


Splitting data table up based on separate file of gene names

There are plenty of efficient and smarter solutions, but if you have so small tables, you can just use the very basics of R to do what you need:
# you need to adjust the following two lines to fit to your file formating and file names
genes_filt_table <- read.table(..., col.names = c('geneid'), stringsAsFactors = F)
norm_counts <- read.table(..., col.names = c('geneid', 'norm_counts'), sep = ',', stringsAsFactors = F)

genes_to_filt <- norm_counts$geneid %in% genes_filt_table$geneid
filt_genes <- norm_counts[genes_to_filt, ]
non_filt_genes <- norm_counts[!genes_to_filt, ]

write.table(filt_genes, 'filtered_genes.csv', quote = F, row.names = F, sep = ',')
write.table(non_filt_genes, 'non-filtered_genes.csv', quote = F, row.names = F, sep = ',')

This being said, there are like milions of more efficient solutions out there. I provided one that is using as basic R functions as possible, as it seemed that you would like to see a simple R solution.


The codes are usually available at DCC's GitHub.
You can always email them if you have more questions.


How to determine the primary Uniprot accession number from a set of accession numbers?

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

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

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

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

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

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

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

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

Both reviewed and unreviewed entries can have primary accession numbers.

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

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


How to extract gene expression tables from this GEO dataset?

Assuming that you have downloaded processed file, according to your link, it is either plain text (compressed) or GTF. If you would have plain text you probably would have posted a snapshot or sample of it, so I assume you have the GTF. According to format description for GTF-files:

Format consists of one line per feature, each containing 9 columns of data, plus optional track definition lines.
The 9 columns (or fields) are tab-separated and all fields contain a value (missing/NA-values are usually denoted with ".").
Fields are:

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

You can find more information on it here and here.

IF you don´t have the file and just have the query output, you probably will find the answer here, where extraction of expression data has been discussed in more detail.


What is the standard way to measure contig sequence lengths in a BAM?

Standard approaches would be either samtools idxstats or samtools view -H, both of which will produce the exact same results. In fact, all methods will produce the exact same result, since contig length is set by the reference you align against.


How to calculate Centimorgan lengths from run lengths in one-to-one dna profile comparison

Try having a look at this link:

http://bochet.gcc.biostat.washington.edu/beagle/genetic_maps/plink.GRCh37.map.zip

If you unzip the file, there is a file for each chromosome (1:22) and then 3 for the x chromosome. The first column of each file details which chromosome the file corresponds to, the 3rd column gives you the cumulative genetic distance in cM and the 4th column gives you the genomic coordinate (position) in bp.

You can then calculate the distance in cM between two positions from this data by subtracting the cM positions between the two positions you are interested in.


How can I ask snakemake to produce a dag where each node represents a rule?

I think you are looking for --rulegraph instead of --dag

  --rulegraph           Do not execute anything and print the dependency graph
                        of rules in the dot language. This will be less
                        crowded than above DAG of jobs, but also show less
                        information. Note that each rule is displayed once,
                        hence the displayed graph will be cyclic if a rule
                        appears in several steps of the workflow. Use this if
                        above option leads to a DAG that is too large.


Can I index a compressed FASTA file using STAR?

I can think of two possible workarounds. The simplest would be to try using zcat instead of gunzip -c, in case that works. Assuming it doesn't, you can use process substitution, if your shell supports it. Something like this:

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles <(zcat ref.fa.gz) --sjdbGTFfile ref.gtf --sjdbOverhang 100

The <() creates a file descriptor and passes that to the command. It is very useful for passing streams to commands that expect files.

If your shell doesn't support this construct, assuming you are on a system that supports named pipes (Linux does), you can do the same thing a bit more manually:

Create a named pipe (FIFO)

mkfifo foo

Set it to read the uncompressed genome

zcat ref.fa.gz > foo &

The & is important since the command above will hang until something starts reading from foo. Alternatively, you can just run this in one terminal and the next command in another.
Read the input from the pipe

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles foo --sjdbGTFfile ref.gtf --sjdbOverhang 100


Hardware Requirements (specs) for Bioinformatics-dedicated desktop

The RAM requirements for some bioinformatics analyses like assembly can be quite high (in the hundreds of Gigabytes). My recommendation is to get a fast laptop. Something with an i7 quad-core processor, 16 GB of RAM, and 1TB of storage should do. Do your analyses on this laptop. If you run out of RAM or storage, then use a cloud service to get access to bigger machines or buy an external hard drive. Once you have a feeling for what your actual needs are, start looking at a dedicated rack server solution in your local institution's server farm or again, in the cloud.

You don't mention any analyses that clearly benefit from a GPU, but as you get more familiar with your analyses, perhaps that will become useful as well.


OK, I got an answer here. Use:

frank off

to get:


Retrieval of genomic position by using biomaRt package

BiomaRt cannot do that I'm afraid. There is an Ensembl REST API endpoint that will get you the genomic location of protein coordinates. It needs an Ensembl peptide ID an input though, so you could use the xref endpoint first to get that. There's a bit of sample code in R on those pages, you could use that to script together if you've got a list of these.


So, columns 7 and 8 should be the genotype of variant 1 (each column is one of the alleles). You have them as 0 1, which means absent and major allele.

I think you're just missing a column. From the doc:

Column 6 is Phenotype value ('1' = control, '2' = case, '-9'/'0'/non-numeric = missing data if case/control)

Ideally you would know which SNPs are heterozygous based on your SNP calling and have such SNPs as 1 2, but if you have a homozygous organism or if you're working with data that's been flattened you can simply generate the .ped file with the same number twice for each SNP.


Best practice to perform a functional analysis on a set of VCFs from a cohort of patients

As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.

This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.

Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.


Why do these NCBI representative genomes for ape species have no Y chromosome?

tl;dr: technical difficulties, or sex

For the gorilla, that genome is from a female gorilla (Kamilah the gorilla) so she doesn't have a Y chromosome.

For the chimp link, there is a visualization of the chromosomes that you can click on that takes you to the genome browser for the Y chromosome (see image). Notably, there is a more recent genome build that does have a Y chromosome sequence assembly that is directly linked, they just haven't updated the main reference sequence yet.

To directly address the question, the Y chromosome is very hard to assemble. The easy part of it to assemble is already on the X chromosome anyways (the pseudoautosomal region), and the hard part of it is often a mess of repeats, especially in mammals. Many groups will decide to just assemble the X as well as they can and call it good, due to the difficulties associated with assembling the Y. If they are working with a male individual the Y sequence is still there, it just may not be anchored to a chromosome sequence, so it's only in the unscaffolded/unanchored assembly ("master WGS" I believe, usually). It is also probably highly fragmented (it isn't a whole chromosome of sequence stitched together, but rather a lot of little pieces of sequence).

In some cases they will focus all sequencing efforts on female individuals specifically to make their assembly job easier, because in females the X is at the same copy number as the autosomes and doesn't have interference from homologous regions on the Y.


How to find genes whose expression vary the most in time-series?

Maybe my thoughts are too simple, but wouldn't the computation of the variance per gene suffice for this? You really do not need a statistical test, since they are not talking about DEGs, just high-variant genes.


Plotting coverage of annotation over collection of region

If one assumes that repeat elements of a given type (e.g., LINEs) don't overlap each other, then the following will work:

Split your BED file by repeat element, such that you have a LINE.bed, SINE.bed, etc.
Convert those to bedGraph (e.g., awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"1.0"}' LINE.bed > LINE.bedGraph).
Use UCSC tools to convert those bedGraph files to BigWig.
Install deepTools and run computeMatrix reference-point -b 2500 -a 20 -S LINE.bigWig SINE.bigWig LTR.bigWig -R Regions_Of_Interest.BED -o foo.mat.gz
Make a profile plot with plotProfile (plotProfile -m foo.mat.gz -o foo.png --perGroup)

If you do have overlapping regions you'll need to first make a disjoint set of intervals (you can use bedops for this), find the coverage of them (e.g., bedtools coverage) and then continue on with that.


Can the canu assembler output a fastq file of the final assembly just like HGAP4?

I've created an issue on the Canu github repository for this. I'm not aware of any existing functionality to output FASTQ files, but think that this would be a useful feature to have.

It would be possible to create something like this by aligning the original reads to the genome, generating a pileup from that alignment, then determining the sequence likelihood based on that pileup. Unfortunately, I don't have any scripts that I've written to do that... yet.


Is it conflicting: significant difference according to Wald test but not in PERMANOVA and ANOSIM to compare species abundances of metagenomic samples?

Yes, it does make sense to perform differential abundance testing after one doesn't find compositional differences with PERMANOVA / ANOSIM.

It seems to me you are drawing an (implicit) analogy to ANOVA and post-hoc tests like the Tukey Honest Significant Difference test. Intuitively, the analogy makes sense, however, these situations are quite different. PERMANOVA / ANOSIM and DESeq2 are operating on different aspects of the data, with different transformations and under different assumptions.

PERMANOVA and ANOSIM do not operate on the raw counts directly, rather, they test differences between groups by calculating some distance or dissimilarity between the samples. One can choose between several distance / dissimilarity measures, and the test result may be impacted by the chosen measure.

DESeq2 (and other packages), on the other hand, test for species abundances under a parametric model, operating on the raw counts directly.

Thus (sadly) one is not a natural, optional follow-up of the other, they are complementary and will tell different stories about the data.


How to find the sum of all sequences in a multifasta file

grep -v ">" multi.fasta | tr -d '\n' | wc -c

grep -v ">" find all non header lines (i.e. without ">")
tr -d '\n' delete all newlines
wc -c count all characters


Those odds ratio aren't calculated from the raw allele frequencies. They're calculated from coefficients estimated from complex linear models with a bunch of covariates that are run on subsets of the data then meta-analysed together.

It is a bit surprising that the odds ratio is less than one when the frequency is higher in cases, but the p-value indicates that the OR is not significantly different from 1 so it may be down to statistical noise. It's possible that within populations one allele is more common in cases than in controls but that cases are primarily sampled from populations where that allele is uncommon, causing it to be more common in controls overall.  See here. Alternatively, 1/0.98 is pretty close to the value you calculated, so it's possible that the reference allele got switched around at some point.


I've just come across a similar issue. It seems the mouse and human mitochondrial genes use mixed capitalizations. See the following PR for details:

https://github.com/Hoohm/dropSeqPipe/pull/43

Bottom line is, if you ignore the case in the grep pattern it might work (assuming you use human/mouse):

mito.genes <- grep("^mt-", rownames(shox2crepos_mtmg@data), value=TRUE, ignore.case=TRUE)


Why does SBOL have both Component and ComponentDefinition?

The key to understanding this is to think about a system with more than one copy of a genetic construct.

For example, let's say you have a construct that uses two copies of the pBAD promoter, each regulating a different coding sequence. Both of these promoters are pBAD promoters, with the same sequence and type information, so they should share the same ComponentDefinition.

But there are two copies. In the context of the larger construct, we need to be able to talk about each copy separately. For example, we might need to say that the first comes after a MoClo cut site, while the second comes after a terminator. For that, we need the Component object, to have a way to say "this particular copy of pBAD", as opposed to "pBAD in general", which is what the ComponentDefinition tells us about.


How do I interpret the function O(f^-2 log n) in Shekhar et al. (2017)?

Big O notation means "of the order of". Thus $\mathcal{O}(f^{-2} \log{n})$ means that the number of trees required goes with $\log{n}$ for number of species (i.e. if you square the number of genes, you should double the number of trees) and with the inverse square of the shortest branch length (i.e. if you double the shortest branch length, you quarter the number of trees needed).

Big O notation generally allows to to compare two input sizes, but cannot really tell you anything about a single point. I.e. it is hard (impossible?) to tell from this function whether you have sufficient loci.


ENA minion 2D data doesn't actually have 2D data

Probably not. ENA doesn't have any ability to specifically flag a file as 2D, and the nanopore file formats have changed a lot, so you pretty much need to rely on the submitter to separate out the reads in a meaningful way.

Also, there wasn't really any consistent convention for what ONT used in their read headers. Looking at some of my old files they had a plain sequence ID (i.e. without anything suggesting temp/comp) in different template/complement folders in the FAST5 files for the components of a 2D read and the consensus sequence, but later changed it to _2D_template and _2D_complement as in the example you have (possibly with an intermediate _1D_template and _1D_complement format). Parsing/understanding those FAST5 files can be difficult, which is why poretools was created in the first place.

It's still a bit of a mess, with the 1D$^2$ reads ending up in a Frankenstein-like state that's partially 1D and partially 2D.


Filter Trinity transcriptome based on RNASeq reads

I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:

Run BUSCO in short/transcript mode on the Trinity-generated sequences
Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
Sub-select the transcripts based on this threshold, using the Salmon results table

Step 1 -- BUSCO

python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10

Step 2 -- Salmon

~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A

Step 3 -- Results merge (using R)

#!/usr/bin/Rscript
busco.table <-
    read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
               header=FALSE, stringsAsFactors = FALSE,
               col.names=c("Busco id","Status","Sequence","Score","Length"));
busco.table <- subset(busco.table, Status != "Missing");

count.table <- read.delim("quant/reads-all_quant_ISR/quant.sf",
                          stringsAsFactors=FALSE);

busco.expr.df <- merge(busco.table, count.table,
                       by.x="Sequence", by.y="Name", all.x=TRUE);

Step 4 -- Count thresholding (visual exploration)

pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
options(scipen=10);
busco.numReads <- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                              function(x){
                                  exp(mean(log(x[x>0])))}), decreasing=TRUE);
par(mfrow=c(1,2));
plot(x=seq(0,1,length.out=length(busco.numReads)),
     busco.numReads, log="y", ylab="Number of Mapped Reads",
     yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
     xleft=0.9, xright=par("usr")[2],
     col="#80808020", border=NA);
plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
     busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
invisible(dev.off());

Looks like a value of 50 will be suitable.

Step 5 -- Sub-select the transcripts

cat(file="HighCount_Transcripts_Num50.txt",
    subset(count.table, NumReads >= 50)$Name, sep="\n");

Then filter using a suitable FASTA selection utility:

pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta


Merging regions according to their identifier

If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g.

awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \
| mergeBed \
| sed 's/_/\t/' \
| awk -v OFS='\t' '{print $1, $3, $4, $2}'

(Assuming "_" is not part of chromosome or identifier name, adjust accordingly if so).


How can I find out if my gene of interest is duplicated?

You might want to have a look at What is the difference between orthologs, paralogs and homologs?. Paralogs are duplicates by definition, and whether a gene is a pseudogene or still active has no bearing on whether it comes from a duplication.

In any case, to figure out what's going in in your target genome, you should use the target's genes. So, first get the best hit for each of your query genes in the target genome. Then, to investigate its evolutionary history, use that sequence (the one from the target genome) as  query against the target genome. That will let you assess how similar the target's hits are to one another and that's how you can attempt to figure out their evolutionary history.

Note that 20% identity is very low for closely related species. While there is no golden threshold, I would expect recent duplication events (and in closely related species, a duplication event that is in one genome and not the other will be recent) to produce much more similar sequences. In any case, you can't do this based on sequence identity alone.

Are the protein domains of the gene products conserved? Is the exonic structure of the gene conserved? You should also try taking the putative duplicate of geneA from genome2 and blast it against genome1 (your original source genome). Is the best hit the geneA from genome1, the query you used to find the putative duplicate in genome2? Or is there a better hit for it against another gene?

These are all questions to investigate and take into account when reconstructing the evolutionary history. You can't just look at identity and use some threshold to decide this sort of thing.

If you blast the potential dupicate from your target genome against your original source genome, is the best hit the duplicate?


RNA seq .counts.txt to bigwig conversion

No, count-level data doesn't have the base-level resolution that is required in order to be usefully visualised on a genome browser. You need a BAM file (or raw FASTQ files that can be used to generate a BAM file) in order to do that.
"No" isn't completely correct: it'd be possible to take gene counts and information about gene annotation boundaries to shoehorn them into a BedGraph / BigWig file, but whatever comes out the other end of that isn't going to be a good representation of reality.


I use the DropletUtils package from Bioconductor. Getting counts begins like this:

library(DropletUtils)
filePath <- "path_to_data/raw_feature_bc_matrix/"
singleCellExperiment <- read10xCounts(filePath, col.names = TRUE)


RNA-Seq data transformation prior to sample correlation analysis

You don't absolutely need to perform any further steps before calculating the correlation. However you would likely benefit from the following sort of workflow:

Transform your data with vst() or robust log (see the DESeq2 vignette) so genes are on a more or less similar scale.
Select the top 300-500 most variable genes
Use this subset with to compute sample-wise correlation.

This also what follows what you would find in the DESeq2 user guide, which I strongly encourage you to read.


I m not sure about htseq-counts but I would suggest go for featureCounts, I will update my answer and share the script/command. Give it a try

the command i use

featureCounts -T 50 -C -p -t exon -g gene_id -a /run/media/data2/gencode.v21.annotation.gtf  -o PRIMARY_CELL.txt $(ls *.bam | sort )


Blast hits disappearing after changing -evalue

I am the author of the question / answer here Why does a very strong BLAST hit get lost when I change num_alignments, num_descriptions or max_target_seqs parameter?

  To my knowledge changing the evalue threshold should no change the
  values themselves, it should just change the reporting.

That is sadly not the case, though thats what I would say is typically expected.

Yours was the exact situation that lead me to question the implementation of NCBI blast. I got results with an evalue above a certain value with a relaxed threshold, but with a stricter threshold (still lower than the value of the results I got) those results disappeared.

From what I found all parameters imposing some threshold show this behavior and are applied in an early step of the program leading to prematurely removed results that later on might score higher than those thresholds. (Say your hit first has an evalue of 5: threshold 10 = ok and kept, threshold 1e-5 = gets removed, but that same hit gets improved further along in the program and then gets an evalue of 1e-6. Now this is only shown when you use an evalue threshold of 5 or higher, though 1e-6 is clearly also a result you would want to see when using a threshold of 3, 1, 0.01, ... etc.)

I was relatively shocked by this and the fact that its basically nowhere officially documented and generally a not known fact what influence these options have.

Sadly for now the stance of the developers is basically all is working as intended and the only feedback I got amounts to

  The safest approach is to relax (increase) the threshold and deal with a larger result set.

aka Suck it up and deal with it.


Ultimately my colleague helped me to solve the issue by following the steps:

Created environment: conda create --name myenv
Activated it: source activate myenv
Edited .condarc to set the priority of the channels to be the following (top to bottom): conda-forge, bioconda, r, default
Checked that the channels are set in the correct order: conda config --get channels
conda install r=3.5.1
conda install r-essentials
conda install rstudio
Installation of all of the bioconductor packages through bioconda and never install anything from inside R or Rstudio (I guess that is how it was broken)

Then launching rstudio from within the environment. Maybe it is not a kind of quick fix and it could have been solved easier, but I was not able to.


How to extract / convert gff3 CDS sequences to multifasta

I like bedtools getfasta.  My typical option set is bedtools getfasta -fi <reference> -bed <gff_file> -name -s.  Be aware of the -s to make sure you are pulling the correct strand.  I like bedtools because it is a versatile tool overall for handling bed, gff and vcf file manipulations.

# bedtools getfasta

Tool:    bedtools getfasta (aka fastaFromBed)
Version: v2.26.0
Summary: Extract DNA sequences from a fasta file based on feature coordinates.

Usage:   bedtools getfasta [OPTIONS] -fi <fasta> -bed <bed/gff/vcf>

Options:
    -fi Input FASTA file
    -bed    BED/GFF/VCF file of ranges to extract from -fi
    -name   Use the name field for the FASTA header
    -split  given BED12 fmt., extract and concatenate the sequencesfrom the BED "blocks" (e.g., exons)
    -tab    Write output in TAB delimited format.
        - Default is FASTA format.

    -s  Force strandedness. If the feature occupies the antisense,
        strand, the sequence will be reverse complemented.
        - By default, strand information is ignored.

    -fullHeader Use full fasta header.
        - By default, only the word before the first space or tab is used.


The plot is probably being drawn to a device that you cannot view on a cluster. If you're using the plotNetworkHeatmap function, say, try something like this:
png("plot.png", ...) #specify height, width, etc
plotNetworkHeatmap(....) #call your function
dev.off()

Now, scp the plot.png file over to your local machine to view the plot.


Why are there more barcodes than GEMs in 10X chromium data?

Your initial guess is almost certainly correct. I don't know about the linked read libraries, but in the 10X single cell sequencing protocol, separating real barcodes from noise barcodes is an important and sensitive step in the analysis pipeline.

The current analysis pipeline is to look at a plot, like the one you produced and find the "knees". The plot you produced has two knees, on at around 10 and one at around 200. I'd guess that the ones below 10 are what we call "error barcodes" generated by sequencing is PCR errors from real barcodes, and the ones between 10 and 200 are what we'd call "background" barcodes - in scRNAseq they represent things like GEMs that didn't encapsulate a cell, but did encapsulate some non-cellular DNA or RNA.

As there are more than 1,000,000 possible barcodes, not all of them are used.  I believe 10X knows which sequences are used and which are not (ie the barcodes are not random), so I'm surprised that these sequences make it out of the other end of the analysis.


Blueprint does not have technical replicates, you will find that very few sequencing based experiments ever have technical replicates. As biological replicates: each individual is regarded as a replicate.

This makes sense for the sort of experimental design which blueprint was designed to answer questions such as "Do individuals from group A have higher or lower H3K4me3 than individuals from groups B" or "Does cell type A have higher accessibility than cell type B". In each of these cases the within group variability to be considered is encapsulated in the difference between individuals.

Overall technical variation is assumed to be similar in all samples or at least not correlated with the contrast of interest (this may well not be true in the second example design I gave) it is safe to directly measure it, as long we know its distribution.

As an aside, one of the reasons that ENCODE got away with only two reps is that their reps are effectively technical (repeated preparations from the same cell lines) and thus there is next to no biological variation (which is what our algorithms are set up to find).


Finding orthologues using BLAST on the NCBI database

BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species.

Since you are new to bioinformatics I recommend you search an ortholog database that already exists.

Here are some databases:

orthoDB
OMA
EggNog

All of these should allow you to search either with text or a sequence to identify the orthologous groups that your input belongs to at different taxonomic levels. Sequence searches either follow pairwise alignment (e.g. BLAST, DIAMOND, RAPSEARCH) or profile HMM type (e.g. HMMER3) approaches.

A more complete list of ortholog databases is available from "Quest for orthologs".


How to scale the size of heat map and row names font size?

heatmap.2 is very configurable, and has options to adjust the things you want to fix:

cexRow: changes the size of the row label font.

keysize: numeric value indicating the size of the key.

The size of the key is also affected by the layout of the plot. heatmap.2 splits your plotting device into 4 panes (see the picture below), and you can control the size of the key partly by controlling the size and layout of these 4 panes.

The layout and scale of these panes can be controlled using the lmat, lwid and lhei arguments. For example, in your case, because you're not displaying the column dendrogram in your plot, you could set the height of the first row of the layout to be smaller, which would reduce the size of the colour key: lhei=c(2, 10) (this is just an example, you'd need to experiment to find values which worked well for your specific heatmap).


Why is 42 the seed used in RunPCA() in Seurat?

Given that it is virtually impossible for a human to predict what random numbers will be generated given a certain seed, and what effects they will have for a given application, the choice of the seed can be considered arbitrary.

What matters is that fixing the seed should fix the behaviour of the application. What seed has been chosen does not really matter, as long as it is documented.

The choice of 42 could therefore well come from one of the reasons you propose. Only the authors of the programs could confirm.


How to find possible viral subsequences in a Eukaryote genome?

Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.


probeset to probeset mappings between Affymetrix arrays

If your question is: can probeset IDs from different platforms be mapped to one another in a similar way as mapping probesets to genes, then the answer is: Yes. BioMart allows you to map almost anything that has an ID to anything else that has an ID.

You can use BioMart either via the web interface or programatically. A brief guide to using the web interface, for mapping HG U133 Plus 2 to HuGene 1.0:

Go to the start page
Select H. sapiens Ensembl genes for your database; Human genes (GRCh38.p10) for your dataset
Click Filters in the left-hand column
Expand "REGION", scroll down to GENE and select "Input microarray probes/probesets ID list [Max 500 advised]"
Select AFFY HG U133 PLUS 2 probe ID(s) and either copy/paste or upload a list, one per line
Click Attributes in the left-hand column
Scroll through, uncheck what you don't want to see and choose what you do, for example Gene Stable ID, AFFY HuGene 1 0 st v1 probe and AFFY HG U133 Plus 2 probe
Finally click "Results" in the menu at the top of the left-hand column

And you should see a result like this (you'll need to click the "Results" button).

You should not expect that there be a one-to-one mapping, for two reasons:

Genes have multiple transcripts and probesets are mapped to each transcript
Some transcripts have more than one probeset: in this case, the HuGene IDs 8002301 and 8002303 map to transcripts for this gene

Finally: here's a programmatic example using R/BioMart:

library(biomaRt)
mart.hs <- useMart("ensembl", "hsapiens_gene_ensembl")
results <- getBM(attributes = c("ensembl_gene_id", "affy_hugene_1_0_st_v1", "affy_hg_u133_plus_2"),
                 filters = "affy_hg_u133_plus_2",
                 values = c("210519_s_at"),
                 mart = mart.hs)

results
  ensembl_gene_id affy_hugene_1_0_st_v1 affy_hg_u133_plus_2
1 ENSG00000181019               8002301         210519_s_at
2 ENSG00000181019               8002303         210519_s_at


low-memory high-speed replacement for Picard MarkDuplicates

As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.


Best way to detect long insertions in bisulfite sequencing data?

You'd be best off by starting with -rfg and -rdg as is and reverting bismark's change of --score-min back to the default for bowtie2. That alone will allow for much longer indels. If that still doesn't suffice, then I'd play around more with --score-min before messing with the gap open/extend penalties. If you do need to play with those, then increase the gap open penalty and decrease the gap extension to 1. Do check that the resulting alignments aren't nonsense though!


How to allow "half windows" in biopython?

You can choose between a window_size of 20 for sequences longer than 20 AAs, and the length of the sequence for those shorter::

seq = "CYGTGAVYLLGLGIGGFS"
window_length = min(len(seq), 20)

full_sequence_analysis = ProteinAnalysis(seq)
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)


Using pythoncyc to get ec-code of a reaction

BioGeek is correct, however I will post a full answer.

import pythoncyc
# To get ec-code
db = pythoncyc.select_organism('sco')
reaction = db["RXN1A0-6128"]
ec_code = reaction["ec_number"]

# To get genes
pgdb = pythoncyc.PGDB("sco")
genes = pgdb.genes_of_reaction("RXN1A0-6128")


Performing PCA for the samples and for the genes

First of all, I am not a statistician but have been performing RNA-seq analysis for a while, here is my take:

I think performing PCA on the samples makes sense in terms of mathematics. What I mean is you can compute eigenvectors and subsequently principal components either on your genes or on your samples.

I don't think performing PCA on the samples makes sense in terms of statistics in RNA-seq context (and probably in any other context): Your genes (features) describe your samples and not the other way around. So the computed PCs should describe/explain your samples and not the other way around. You would not perform PCA on your samples (individual cells in scRNA-seq) even in the case of scRNA-seq, where number of samples sometimes exceeds the number of genes, there PCA is performed on genes and the resulting few PCs are used for downstream analysis, usually clustering and visualization, both used to infer about samples (cells).


How to generate a weblogo image that takes into account sequence frequencies?

Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.

Generate a simple mapping file (counts) with the number of times each sequence was found:

seq1 3
seq2 10
seq3 2

Then, given an aligned protein fasta file (foo.aln) like this:

>seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

You can simply repeat the sequences as needed (see my answer here for the FastaToTbl and TblToFasta commands):

$ awk '
        (NR==FNR){
            count[$1]=$2;
            next
        }
        {
            for(k=1;k<=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts <(FastaToTbl foo.aln) | TblToFasta
>seq1.1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.2
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.3
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2.1
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.3
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.4
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.5
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.6
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.7
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.8
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.9
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.10
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3.1
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
>seq3.2
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

Note that the above assumes your shell supports the <(command) syntax for process substitution. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run FastaToTbl foo.aln > foo.tbl first, and then use foo.tbl as input.


Finding differentially expressed genes between two samples

If I understand your question correctly, you can simply use SetIdent() to change the "default identity" to your samples and then use FindMarker() with the ident.1 = "Double-KO" and ident.2 = "Shox2-KO".

Alternatively, you can "collapse" or "summarize" your single cell data into a pseudo-bulk and then use one of the bulk RNA-seq packages.


What are the labels in my rstats distance matrix?

Disclaimer: I don't have adagenet and hierfstat installed, I just looked at the source code. It seems like pairwise.fst (the deprecated code from the adagenet package, which you seem to be using) already returns a matrix. Set res.type="matrix" and try calling rownames(f) and colnames(f) afterwards.


Installing RnBeads via bioconductor - .onLoad failed in loadNamespace(), call: NULL

I have reported the error on support.bioconductor.org but we were not able to solve the problem. It seems that some dependencies are making problems, but as reinstalling my home R folder did not resolve the issue, I have come up with the following workaround:
I downloaded the latest R version from my favorite CRAN mirror, opened an R session with pathToNewR/bin/R from my terminal and installed RnBeads there. This works fine.
Note: I had to make sure that pathToNewR/library is the only entry in .libPaths(), as per default, my home R folder was also contained in .libPaths().
As a consequence, I am independent from the (maybe broken) R version and packages on my server and can finally work with RnBeads :)


Question about the dots on Quartile groups in boxplot

You add the points with geom_point(). Just remove it and you will get your "empty" boxplot.

q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))

Unfortunately I couldn't use stat_compare_means(method = "t.test") and     stat_n_text(), but this two just add the labels.


If your goal is only to account for degree distribution (i.e. Matteo Ferla's comment), a simulation approach you could use is:

measure the degree (number of links) of each protein in your PPI
find a large set of proteins that are equal or similar in degree to each of your proteins of interest (you may have to use a range of degrees to get a large enough set for comparison, ideally >=1000). We'll call these "comparable proteins" for each of your proteins of interest.
For each protein of interest, make all of your measurements also for its set of comparable proteins
Compare your measurement (shortest distance?) for the protein of interest to that for its comparable proteins. You could compute Z-scores from the parameters of the comparable protein distribution, or just empirical Monte Carlo p-values from where it falls in the distribution.

This will give you at least a heuristic null distribution for each protein. It should be relatively easy to code up in whatever your favorite language is.

Hope that helps-


How to predict stop codons in Illumina reads?

With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.


What unit do I get on the y-axis of a metagene profile plot?

Feel free to @ me in deepTools questions, since I'm the primary developer.

For a given bin, the count assigned to it is the number of reads overlapping it, regardless of whether they overlap by 1 or 10 bases. So a read overlapping only partially and one overlapping completely are treated the same.

Since your bigWig file is in units of "alignments" (i.e., it's not 1x normalized), the resulting profile will also be in units of "alignments" (i.e., profiles and heatmaps are in whatever units the input files are in).

Upstream/downstream regions and unscaled regions are also 10 base bins. Note that these are then the average of the per-base value, since the bins here may not perfectly correspond to the bins in the bigWig files. The line in the profile plot is indeed the average (by default, you can choose median, max, min, etc.) of the underlying regions for each bin.

Regarding the scaled section in the middle, the the number of genomic bases per bin is changed such that the region will have "length"/(regionBodyLength/binSize) bases each. As above, the per-base value is then averaged (or whatever you specify) to derive the per-bin value. The length here is decreased if you have unscaled regions, since otherwise bases would get counted twice.


Which tools for differential expression analysis in scRNA-Seq?

According to Jaakkola et al. (ref. 1):

  Based on our comparisons, DESeq and Limma without any modifications are not suitable for scRNA-seq data analysis, and yet, they have performed well in the context of bulk RNA-seq data.

They also provided a summary of the comparison between 5 tools (2 designed for Single-Cell, 2 designed for bulk, 1 general):

In my case I ended up using MAST (ref. 2, github).

Jaakkola, M. K., Seyednasrollah, F., Mehmood, A. & Elo, L. L. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics 18, 735–743 (2017).
Finak, G. et al. MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 1183 (2015).


size of the pathways for analysis and filtration

The problem here is that any one GO term with a small number of members is unlikely to include any differentially expressed genes by chance, but if you test a large number of GO terms, some of them definitley will. The solution is to correct the p-values for multiple comparisons. If you do this, then none of those GO terms with a single DE member will be significant. The larger GO terms with multiple DE members in the chart you show probably won't survive correction either. Hopefully you have some that are more significant.


How to search a specific sequence in BAM files for 10X experiment

The bam file has bam tags which say what reads belong to what cells.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/output/bam


How to measure or assign hydrophobicity score values to individual amino acids of a PDB structure?

No, this doesn't makes sense (I think).

The hydrophobicity of an AA is the hydrophobicity of the AA, in whatever context, just like it's molecular weight. It's weight is it's weight, in a structure or in isolation.

I know this seems physically unrealistic, but so is an abstract measure of 'hydrophobicity', it's just the best we can do based on partition coefficients.


CRISPR/Cas9 screen analysis with Mageck: paired-end sequencing

The other answer is incorrect. The original poster is asking about paired-end sequencing data, not paired treatment/control experimental design.
Here is how you analyze paired-end sequencing data (this was an answer I gave on Biostars last year):
https://www.biostars.org/p/303656/


Anyone can declare anything to be a "gold-standard", the term is meaningless. All if actually means is, "we think this is the test that everything else should be benchmarked against".


The solution was to use a different library as I expected. It should be library(TxDb.Hsapiens.UCSC.hg19.ensGene). And also, as specified in the question, to remove everything after the dot in the Name column of the quant.sf files.


Retrieving gene matrix from internet in GSEA desktop

Have you checked firewall settings on your laptop?

See Firewall / FTP connection issues at:

GSEA wiki - known issues


To expand on the excellent answer by haci:
BLAST reports (a) the summary of alignments and (b) individual alignments, see below. Note that your link to the results expired (as NCBI BLAST results links tend to do). It is considered best practice to post the text version of the relevant portion of the results in the question itself.
The summary provides an overview of all the query sequence individual hits in the subject sequence. Because the individual hits (alignments, or ranges) together cover 100% of the query and each of these has the max identity = 100%, that is what it reports in the summary. It is not very intuitive for many users.
You should always look at the individual hits in detail. The above 100% coverage in the summary consists of ranges 1 and 14, shown below. The top hit (range 1) covers only 17 out of 21 nt with 100% identity, and this is the most relevant information from BLAST. The other hit (range 14) covers 10 out of 21 nt with 100% identity. This is the cause of the discrepancy.
These 2 alignments are not close together in the subject sequence. Depending on your goal, choose one or more of the individual alignments and consider alignment length and % identity of those hits individually, rather than summarized together as in the BLAST report summary. For most applications, the top hit (range 1) would be the most relevant alignment.
Program: BLASTN
Subject:Neisseria meningitidis strain M22814 chromosome, complete genome ID: CP031332.1(nucleic acid) Length: 2190201
Query #1: seq1 Query ID: lcl|Query_56011 Length: 21

Sequences producing significant alignments:
                                                                  Max    Total Query   E   Per.
Description                                                       Score  Score cover Value Ident  Accession
Neisseria meningitidis strain M22814 chromosome, complete genome  34.2   57729 100%  0.001 100.00 CP031332.1

Alignments:
>Neisseria meningitidis strain M22814 chromosome, complete genome
Sequence ID: CP031332.1 Length: 2190201
Range 1: 2074883 to 2074899

Score:34.2 bits(17), Expect:0.001,
Identities:17/17(100%),  Gaps:0/17(0%), Strand: Plus/Plus

Query  1        CAATTGCGAAGCCACAT  17
                |||||||||||||||||
Sbjct  2074883  CAATTGCGAAGCCACAT  2074899

...

Range 14: 354292 to 354301

Score:20.3 bits(10), Expect:15,
Identities:10/10(100%),  Gaps:0/10(0%), Strand: Plus/Minus

Query  12      CCACATTATT  21
               ||||||||||
Sbjct  354301  CCACATTATT  354292


Discordance of Gene Ontology data between Mouse Genome Informatics and Ensembl sources

Short answer, we get our annotations from GOA, which has the three listed. We're in touch with GOA to see why they differ from MGI.


This is usually done using AnnotationDbi in combination with the annotation database for human org.Hs.eg.db.
library(org.Hs.eg.db)
genes <- c("ENSG00000102055", "ENSG00000232605", "ENSG00000136828",
           "ENSG00000159197", "ENSG00000201394", "ENSG00000278580",
           "ENSG00000213440", "ENSG00000166816", "ENSG00000275895",
           "ENSG00000105723")
select(org.Hs.eg.db, genes, c("ENSEMBL", "ALIAS"), "ENSEMBL")

This returns:
           ENSEMBL     ALIAS
1  ENSG00000102055       I-4
2  ENSG00000102055  PPP1R2P9
3  ENSG00000102055   PPP1R2C
4  ENSG00000232605      <NA>
5  ENSG00000136828   RALGEF2
6  ENSG00000136828  RALGPS1A
7  ENSG00000136828   RALGPS1
8  ENSG00000159197     ATFB4
9  ENSG00000159197      LQT5
10 ENSG00000159197      LQT6
11 ENSG00000159197     MIRP1
12 ENSG00000159197     KCNE2
13 ENSG00000201394   RN5S140
14 ENSG00000201394 RNA5SP140
15 ENSG00000278580      <NA>
16 ENSG00000213440      <NA>
17 ENSG00000166816     DLACD
18 ENSG00000166816       DLD
19 ENSG00000166816      LDHD
20 ENSG00000275895   U2AF1L5
21 ENSG00000105723     GSK3A

If you want to apply a 1:1 mapping of gene IDs to common names, you can use the "SYMBOL" column. In a data.frame, you can do something like this:
ensembl_to_symbol <- function(e) {
  hits <- select(org.Hs.eg.db, e, "SYMBOL", "ENSEMBL")$SYMBOL
  # Gene IDs without symbols should just stay as IDs
  no_hit <- is.na(hits)
  hits[no_hit] <- e[no_hit]
  hits
}

recode_symbols <- function(s) {
  e <- ensembl_to_symbol(unlist(strsplit(s, split = ",")))
  paste(e, collapse = ",")
}

df <- data.frame(
  genes = c(paste(genes[1:5], collapse = ","),
            paste(genes[6:10], collapse = ",")
  )
)

df$genes <- sapply(df$genes, recode_symbols)

This isn't the most vectorized (fast) solution, but it's probably easier to understand this way.


How to subset a GRanges via an argument passed into a function?

subsetter = function(gr, cname) {
    return(mcols(gr)[[cname]])
}

You can then use things like subsetter(gr, "GC") and subsetter(gr, "score").


Pattern matching and removing from the name

Without regex, split on "_", drop last item, then paste it back again:

sapply(strsplit(x, "_"), function(i) paste(head(i, -1), collapse = "_"))

Or with regex (stealing the code from tools::file_path_sans_ext function which is used to drop file extensions):

sub("([^_]+)\\_[[:alnum:]]+$", "\\1", x)

# [1] "BAF60A_KD_1"     "BAF60A_KD_2"     "BAF60A_RA_1"     "BAF60A_RA_2"     "BAF60A_VD3_1"    "BAF60A_VD3_2"
# [7] "BAF60C_KD_1_UI"  "BAF60C_KD_2_UI"  "BAF60C_KD_RA_1"  "BAF60C_KD_RA_2"  "BAF60C_KD_VD3_1" "BAF60C_KD_VD3_2"
# [13] "VC_1_UI"         "VC_2_UI"         "VC_RA1"          "VC_RA_2"         "VC_VD3_1"        "VC_VD3_2"


I'm not sure I understand the question, but XP_019970915.1 is an 'NCBI Reference Sequence' ID for 'erythrocyte binding antigen-181' (https://www.ncbi.nlm.nih.gov/protein/XP_019970915.1/). If you look that protein up in UniProt you can find other identifiers/info and GO annotations etc (https://www.uniprot.org/uniprot/A0A151LVZ9). Does that answer your question?


Extracting sequences from FASTA beginning with common 5' end

Here is a 'quick' python3 solution, requiring installation of BLASTN and a pip3 install biopython.  The results are output in FASTA format to stdout.

Save the following python3 code in hit_finder.py

#!/usr/bin/env python3

import sys
from Bio import SeqIO
from Bio.Blast.Applications import NcbiblastnCommandline

cline = NcbiblastnCommandline(query=sys.argv[1],
                              subject=sys.argv[2],
                              outfmt="'6 qseqid length qstart'",
                              max_hsps=1,
                              word_size=4)

result = [j[0] for j in [i.split('\t') for i in
                         cline()[0].rstrip().split('\n')]
          if int(j[1]) >= int(sys.argv[3]) and int(j[2]) <= int(sys.argv[4])]

for record in SeqIO.parse(sys.argv[1], 'fasta'):
    if record.id in result:
        print(record.format('fasta'))

Make it executable: chmod +x hit_finder.py.

This script takes four arguments (which in future versions could be enhanced using argparse instead of sys.argv):

arg1 = query_fasta_file (file to search in)
arg2 = subject_fasta_file (database file to contain the sequence to match at 5')
arg3 = min_alignment_length (this means the length of the hit, allows mismatching in string with blastn and a word_size of 4)
arg4 = start_of_hit_at_5prime (where should the alignment hit start?)

Let's say you wanted to search for AATCTGTAGTACGACTCACTATAGC in your query seqeuences.  Save this sequence in subject.fa as:

>subject
AATCTGTAGTACGACTCACTATAGC

Your query sequences are stored in queries.fa as:

>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

Now, to search for hits of at least 23 bp starting no further than 1 base from the 5' end of your sequence, you would run this search using:

./hit_finder.py queries.fa subject.fa 23 1

To search for hits of at least 30 bp no more than 10 bp from the 5' end of your queries, do:

./hit_finder.py queries.fa subject.fa 30 10

To save your results, output stdout to file, e.g.,:

./hit_finder.py queries.fa subject.fa 30 10 > results.fa


How to modify execution orders for Snakemake?

I think you should use the priority directive to give precedence to downstream rule(s).

Here below I give higher priority to step2 so when the first sample has completed step1 snakemake will run step2 on that sample rather than submitting another sample to step1.

SAMPLES = ['S1', 'S2', 'S3']

wildcard_constraints:
    sample= '|'.join([re.escape(x) for x in SAMPLES])

rule all:
    input: expand("{sample}.done", sample = SAMPLES)

rule step1:
    priority: 1
    output: "{sample}.step1.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

rule step2:
    priority: 10
    input: "{sample}.step1.done"
    output: "{sample}.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

(Note that without wildcard_constraints your/this code throws AmbiguousRuleException)


Uploading files from local computer or web in rstudio cloud

This seems to be an encoding problem (certain characters not being recognized), see this, this and also this. Probably the file "/cloud/project/SCG71.loom" is encoded with a format that is not included in your "locale".

Check the output of Sys.getlocale() and see if your input file's encoding is listed. If not, you will need to try something like Sys.setlocale("LC_CTYPE", "en_US.UTF-8") or Sys.setlocale("LC_ALL", "English") (see the third link for details).


Separate multiple sequence into different file, one sequence per file

There's an app for that! The exonerate package (see here) includes a nice fastaexplode command that does exactly what you're looking for and will be much, much faster than anything we give you here since it is written in C and not a scripting language. Just install exonerate (you can do that with sudo apt install exonerate on Debian-based systems) and then:

fastaexplode file.fa

This will create these two files (file.fa is the original input file from your question):

$ ls
file.fa  KQK14546.fa  KQK21959.fa


Marking optical or PCR duplicates with picard vs. samtools flagstat

There are duplicates, in this line:

1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356

According to samtools documentation for flagstat:

  Provides counts for each of 13 categories based primarily on bit flags
  in the FLAG field. Each category in the output is broken down into QC
  pass and QC fail. In the default output format, these are presented as
  "#PASS + #FAIL" followed by a description of the category.

So you can usually take whatever is in the first column


I cannot access the book, it might be an older version or there are some steps in between. The screen shot doesn't help.
If you look at the vignette and the error, what you need to provide is a matrix of characters with defined row names, for example:
set.seed(111)
x = sample(letters,100,replace=TRUE)
phyDat(matrix(x),"USER",levels=unique(x))

Error in phyDat.default(data, levels = levels, return.index = return.index,  :
  data object must contain taxa names

mat = matrix(x,10,10,dimnames=list(paste0("t",1:10),NULL))
phyDat(mat,"USER",levels=unique(x))
10 sequences with 10 character and 10 different site patterns.
The states are n t s o y e c h z q m j d o h n z y x p q g u o a i h y c g t d f y t d f e u e f d i j r z x w u m a j p j g d t q y v o x p m s o z x z a y j z r u s a y k l f g l f r i l r f m t t e f x v f g q j


about the normalization of RNAseq data for calculating distance for unsupervised learning

Typically some sort of variance stabilizing transform is used before clustering. Popular options are regularized log transformation or a vst transform, which are available in DESeq2. Note that these are NOT used for performing differential expression, just things like clustering and PCA. For differential expression one would use TMM, RLE or similar.


You want FASTA files from here.
If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).
See the README for more info.


How to download the whole directory of an ensembl FTP page?

You can use the mirroring function of wget:

wget -m ftp://ftp.ensembl.org/pub/release-91/embl/bos_taurus/


What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?


The problem in your code is that the function atp_binding is not returning any value, but printing it to the stdout.

You should add a return statement with the result you want to save in your output file.

Also, you should change how the resulting value is being saved. As the return value of your function is None, the file handler will write nothing.

I would use something similar to this:

def atp_binding(records,patern):
    s = ""
    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()
            iend=iend+50
            s+="{0}\t{1}\n".format(istart,m.group())
    return s

file = open("output.txt","w")

for seq_record in SeqIO.parse(records, "fasta"):
   file.write(atp_binding(seq_record.records,patern)+"\n")
   file.write('-------------------\n')
file.close()


Adding treatment groups via metadata to Seurat object?

You should have fed "1h" to your function call but you fed 1h instead. The latter is looking for an object with that name and as the error shows, there is no such object.


Running python script from terminal: umap issue

You answered your own question, yet it is not a coincidence that you have to change the name of your script.

When Python executes your custom script, it reads:

  import umap

Python will first search the local sources (the directory of the script) , and afterwards the installed sources.

In this case Python identifies the local source to be umap.py, and imports the (currently executed) script, while foregoing the installed sources.
Since you are calling functions, not defined in your current script, Python will not be able to fully execute your script.

For more information on importing see Module search path at docs.pyhton.org


Access base aligned to particular reference position

You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.

for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

That's not exactly correct syntax, but you get the idea.


Why can I not install snakemake on my SLURM computer even though I can find it in the bioconda channel?

I created a file .condarc in my home directory with the following content:

channels:
  - conda-forge
  - bioconda
  - defaults

and it works now. I don't understand why this is necessary though. I thought the channels in the .condarc file would be used if I did not specify them in my command.


Gene set enrichment analysis on differential phosphorylation sites

Your dataset looks perfect for the SUMSTAT enrichment test.

You need to come up with a statistic representing your gene. The simplest ideas here are the number of sites or, probably better, a proportion of phosphorylated sites.
Now you can compute a statistic for every gene set, for example just sum (SUMSTAT). You can also have average, sum of squares or something else.
Geta a null distribution of your gene set statistics by permutations. You can either randomize statistics for genes, or just randomly assign phosphorylated sites across the genome keeping the total number of sites.

Now you can compute p-value just by comparing your value with the null-distribution.

You have to be cautious of two things:

You are testing multiple gene sets, i.e. you perform a statistical test for every gene set. Therefore correct your significance for multiple testing. I suggest using FDR for this.
Beware of potential biases when doing permutations. The obvious one is gene length, i.e. longer genes have higher chance of getting at least one site. But there can be other ones like GC-content, chromosomes, etc. You can overcome this by using more realistic permutations, or by controlling for the potential correlations. You can get a few ideas on fighting biases from this paper (sorry for the self-advertisement).


bcftools has the view command which can do this. Create a text file (samples) which has the names of the samples you want, one per line, and then:

bcftools view -S samples merged.vcf.gz > new.vcf

Or, if you're running a shell that supports it, you can use input redirection for a single sample name:

bcftools view -S <(echo "sample name") merged.vcf.gz > new.vcf


Correct way to assemble reference transcriptome - what's --samples_file option in Trinity?

--samples_file flag is for DE analysis, you can put info about biological replicates in that file. You say that you want to make a ref transcriptome, so you don't need a samples file for DE analysis.

Why do you down sample? Are you not afraid you will not have enough reads to have a complete transcriptome now?


Excess average estimated identical-by-descent in genotype data

A clear negative association within an inbreeding coefficient suggests convergent evolution driven by shared selective pressure. There is an outside chance of purifying selection, or else balancing selection, but selection is the key word. The phenotype of the locus needs investigating, i.e. outside nearly neutral theory what is constitutes positive selection?
Allopatry caused by genetic isolation and genetic drift is your null hypothesis. However Darwinian positive selection is the alternative hypothesis.
I would be wary against suggesting contamination unless you've a clear suspicion this is occuring. I do accept that in human genetics its often difficult to avoid, particularly in ancient DNA type work. In this case you would need to try and control for it. The argument I would put against contamination is that 'genetic relatedness' is localised within the genome, whereas contamination would be a genome-wide phenomena.
Convergence is common within for example HLA genotypes (T-cell) loci between global populations.
You could perform formal tests for postive selection either within pop gen (the locus will associate with significant homozygosity) or phylogeny (point mutation).


Calculating Unifrac Distance and missing data

To calculate the unifrac distance you need the phylogenetic tree. If you only have it for 500 taxa then you can only calculate it for those taxa.

It is impossible to impute the phylogenetic tree without knowing the sequence on which the taxa is based. But if you know it then you can calculate the phylogenetic tree for the 3k taxa you have.


Hierarchial PCA Clustering with duplicated row names


How can atomic RMS differences / atomic RMS distributions help us understand the structure of a protein?

The "RMS distribution about the mean coordinate positions" could be restated as "certainty of the structure". In other words, the main backbone that defines the large-scale shape in the structure is expected to be accurate to within about 0.41 angstrom, while the atoms more on the periphery (or otherwise not forming the backbone of the structure) are less accurately solved. How important this is will depend on what you want to do with the structure. If you want to get an idea for what it looks like then the accuracy is less important, if you want to do docking then the more accurate the atoms around any sort of pocket/interaction surface are the better.


What is the best distribution to model the FPKM values from normalized RNA-Seq data?

FPKM/TPM values are generally log-normally distributed. Reference : Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels


You can get standardised residuals from chisq.test and convert these to p-values, but it's not recommended. The frequencies of each genotype depend on the counts of all other genotypes, so the p-values aren't independent and they aren't corrected for multiple testing, but if you must:

pnorm(-abs(chisq.test(tmp[,c('A','B')])$stdres)))


First I installed this libraries:

sudo apt-get install libcurl4-openssl-dev libxml2-dev

I then installed RCurl:

BiocManager::install("RCurl")

and lastly:

BiocManager::install("DESeq2", version = "3.8")