RNA seq fasta file annotation from alignment to reference matches

An awk approach (assuming your input "csv" file is space-separated):

$ awk '(NR==FNR){a[">"$1]=$2; next;} (/>/){$1=$1" "a[$1]}1;' foo.csv foo.fa
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA

Explanation

NR==FNR : FNR is the line number of the current input file and NR is the line number of all input. When processing more than one file, the two will only be equal while reading the first file.
(NR==FNR){a[">"$1]=$2; next;} : if we're reading the first file, the csv, save the second field as the value in the associative array a, whose indices are a > followed by the 1st field. I use the > only to simplify matching it later in the fasta file. The next moves on to the next line so the rest of the command isn't run for the 1st file.
(/>/){$1=$1" "a[$1]}: Now we're in the 2nd file. If the line starts with a >, then set the 1st field to be itself, a space and then whatever was stored in the array a for this 1st field.
1; : this is awk shorthand for "print this line".

You could expand the above into:

awk '{
        if(NR==FNR){
            a[">"$1]=$2;
        }
        else{
            if(/>/){
                $1=$1" "a[$1]
            }
        }
        print;
      }' foo.csv foo.fa

If your files are small enough to be sorted easily, you can use the FastaToTbl and TblToFasta scripts I have posted before and then join:

$ join -o 1.1,2.2,1.2 <(FastaToTbl foo.fa | sort ) <(sort  foo.csv ) | TblToFasta
>Contig11 ref|XM_021487230.1|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
GACTCTCTTAAGGTAGCCAAATGCCTCGTCATCTAATTAGTGACGCGCATGAATGGATTA
>Contig12 mamamama|_PREDICTED:_Mizuhopecten_yessoensis_UNC93-like_protein_MFSD11_(LOC110443183),_mRNA
AKHSALKJHDLKAJBDLKAJDBADADKJALDKJBAD


Annotating a .vcf with centimorgan information

I think a useful tool to address your question is QCTOOL v2. Particularly, the option you are looking for is -annotate-genetic-map. This tool has been developed originally for GWAS datasets. It has multiple inputs data format.
If you could provide input data as an example, I will try working on it to provide the necessary steps to be done.
Hope this answers your question.


When plotting networks where only the weight of an edge is known like this one. The relative position between elements of the networks is irrelevant. The only informative part of the graph is (usually) the size of the edges, where is the thicker means more weight. This isn't a PCA, so the distance between different nodes is not relevant.

For those reasons there are several methods to build these networks, some of them give more emphasis to represent better the degree of the network or other characteristics of the networks.

I consider this representations bad for informing about the networks, because the cut off is arbitrary and the distance between nodes is not informative. One could first do a PCA or MDS and then add the edges between the nodes, this way the relative distance between the nodes would be meaningful as well as the edges. Or one could also use a heatmap to represent all the correlations between the nodes.


Separate boxplots for multiple violin plot

The vioplot package comes built in with boxplots. You can download it from CRAN or there are more features (including formula input and separate colours) in the development version on GitHub: https://github.com/TomKellyGenetics/vioplot

devtools::install_github(“TomKellyGenetics/vioplot”)
library(“vioplot”)
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident)

To plot several PCs:

library(“RColorBrewer”)
colours <- brewer.pal(length(unique(pbmc_small@meta.data$ident)), “Set3”)
par(mfrow = c(1, 2))
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident, main = “PC1”, xlab = “Identity”, col = colours)
vioplot(pbmc_small@dr@pca[,2]~pbmc_small@meta.data$ident , main = “PC2”, xlab = “Identity”, col = colours)
par(mfrow = c(1, 1))


Simulating phenotype with the 1000 Genomes Project

I managed to get it done without the gcta64.
I have written a Jupyter Notebook with the steps: https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb


Change sample names while using cummeRbund

Don't know which file was used by cummeRbund for sample names but replacing q1,q2,q3 and q4  with Batis, Blue_silver, Chakwal_50 and Syn_22 for all the files in the directory containing cuffdiff output worked for me.

I used sed language for this as shown below:

sed -i 's/q4/Syn_22/g' *

Or more generally you can use

sed -i 's/sample_name_shown/sample_name_to_change/g' *

After that I deleted the sqlite object created by cummeRbund and again created this object using readCufflinks() in R.
This time sample names are changed as desired:

samples(cuff)
sample_index sample_name sample_name parameter value
1            1       Batis        <NA>      <NA>  <NA>
2            2 Blue_silver        <NA>      <NA>  <NA>
3            3  Chakwal_50        <NA>      <NA>  <NA>
4            4      Syn_22        <NA>      <NA>  <NA>


This python script will remove all stop codons from your fasta file. I called it remove_stops.py

#!/usr/bin/env python3
import sys
from Bio import SeqIO
stop_codons = ["TAG", "TGA", "TAA", "UAG", "UGA", "UAA"]

fasta_it = SeqIO.parse(open(sys.argv[1]), 'fasta')
for fasta in fasta_it:
    name, sequence = fasta.description, str(fasta.seq)
    try:
        assert(len(sequence) % 3 == 0)
    except:
        raise IOError("length of sequence {} was not divisible by three".format(name))
    print(">{}".format(name))
    print("".join([sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3)
                   if sequence[i:i+3] not in stop_codons]))

You can run it like so:

python remove_stops.py stop_codons.fasta > filtered.fasta

This program won't work if there are any fasta records that do not have a length divisible by three. This indicates a sliding translation stop and these sequences should be trimmed individually. Also, all the genes must be in the 5'->3' direction.

It works by looking at each codon [sequence[i:i+3] for i in range(0, len(sequence) - 3 + 1, 3), then just makes sure they are not a stop codon: if sequence[i:i+3] not in stop_codons]


What is the best way to address the question of doublets and multiplets in a single cell RNA seq data set?

Indeed there's no satisfactory doublet/multiplet identification methods as of 2020.
The most reliable way to remove doublets is still manual removal of doublet subclusters guided by marker genes.
Not only trajectory analysis can be impacted by doublets, but clustering and embedding can also be affected for obvious reasons.


Transcription factor targets for Human and Mouse genes

JASPAR is another useful transcription factor profiling database.

  The JASPAR CORE contains a curated, non-redundant set of profiles, derived from published and experimentally defined transcription factor binding sites for eukaryotes. It should be used, when seeking models for specific factors or structural classes, or if experimental evidence is paramount.

http://jaspar.genereg.net/


How to export web NCBI tBLASTn results in table format with many queries?

It sounds like your easiest option is to just install blast on your local machine and set up the databases there. It really will only take a few minutes to blast 6,000 sentences against a yeast genome on a modern laptop. Any solution interfacing with web results or parsing the output will be painfully convoluted. There are executables for blast that will work on Mac, Linux, or Windows.

Some beginning instructions for installing it on your laptop can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download . Make sure to pass the --num_threads option to speed up the search.


How can I reproduce a manual NCBI search with Biopython Entrez module?

This can be done by using the "Search details" as a search term in Entrez.esearch:

import os
OPJ = os.path.join
base_dir = os.getcwd()

from Bio import Entrez
Entrez.email = "my_email_here"
query = "ANOS1[gene name] AND refseq[filter]"
handle = Entrez.esearch(term=query, db="nuccore", retmax=400)
ids = Entrez.read(handle)["IdList"]
with open(OPJ(base_dir, "ANOS1_orthologs_max400_ids.txt"), "w") as ids_file:
    ids_file.write("\n".join(ids))
    ids_file.write("\n")


How to convert a SDF molecular structure file to a pandas dataframe with Python?

It is easily done with RDKit:
from rdkit.Chem import PandasTools

fn = '/mypath/dataset.sdf'

df = PandasTools.LoadSDF(fn, embedProps=True, molColName=None, smilesName='smiles')


What are the differences between "Site" and "CONECT" in a PDB file

how come not every atom in the CONECT section is described in the SITE section?

Because it is not required and not needed.

CONECT is for listing bonds, so that molecular graphics programs know how to draw the structure.

SITE is for annotation of interesting sites. Each listed SITE has a corresponding REMARK 800 that "details its significance". A graphical program may provide a menu with a list of all interesting sites for easy navigation.

Note: PDB is a legacy file format. The primary file format used nowadays by the PDB (organization) is called PDBx/mmCIF. So when looking into such details it is often better to look directly into mmCIF files.


Analyzing proteins based on sequence similarity

I suggest you give DAVID a try. Specifically, their Functional Annotation tool. Just enter your list of protein IDs, and it will return groups of proteins where particular GO functions are overrepresented.

This should at least help you start categorizing your list into functional groups.


Convert sequence in MS doc to fasta or genbank

The CrazyDoc Python package can convert Genbank/Fasta/Snapgene/MSDoc sequences into Biopython records, which can be saved as Genbank or Fasta.
A web interface / demo is provided at EGF CUBA (Collection
of Useful Biological Apps): Convert Sequence Files.
Disclaimer: I'm the current maintainer of CrazyDoc.


There are two potential sources of bias in this design.

We cannot distinguish correlation from causation.

Imagine two cases. In the first, the disease progression is inducing immune response. Later stages will be associated with the higher gene expression levels. In the second scenario, the disease is caused by overexpression of a gene. Later stages will be also associated with the higher expression.

This is typical for observational studies. But I just want to mention that a special care should be taken during interpretation of the results.

If we are not following our individuals, we cannot distinguish correlation and avoidance.

Let's say the disease is lethal in certain cases; the survival is negatively correlated with the disease stage. Now imagine there is a gene, which causes severe symptoms when highly expressed. On the later stages you will only observe those patients in which the gene was not highly expressed. From that you would conclude that gene expression is decreasing with the disease progression. In reality this gene is very important and causal, you just do not have patients which are alive to see this.

This is similar to Wald's studies of aircrafts (Survival bias).

  Researchers from the Center for Naval Analyses had conducted a study
  of the damage done to aircraft that had returned from missions, and
  had recommended that armor be added to the areas that showed the most
  damage.
  Wald proposed that the Navy instead reinforce the areas where the
  returning aircraft were unscathed, since those were the areas that, if
  hit, would cause the plane to be lost.

I think the second point is crucial, and can and will lead to false conclusions. I do not see immediately how to avoid this bias if you study different individuals.

Therefore I suggest that the same individuals are followed for a long time.

There are different approaches you can use later. For example you can have two-stage procedure:

Identify genes which are differentially expressed between healthy (H) and sick (A, B or C).
Build a linear model of disease stage stage ~ gene1 + gene2 + ..., using genes identified at step 1.
Similarly build a linear model of survival as a function gene expression.
It is possible you can use logistic regression to infer probability of transition from stage B to stage C.

This is not the only possible approach of course. As you suggested, Markov model is also applicable.

While following the same individuals it is possible to use continuous-time Markov model for such a study. In this case discrete or continuous measurements are recorded at certain time points, and the model parameters are inferred using maximum likelihood.

I'm not an expert in this field, but I think this paper describing the msm package for R will be useful. It also supports hidden Markov models.


You can express the annotations of the genome with a dot plot, x and y represent the two annotations along the genome. The you can summary the result with Spearman correlation coefficient(s) based on whole set of annotations or particular gene families.


Fast way to count number of reads and number of bases in a fastq file?

It's difficult to get this to go massively quicker I think - as with this question working with large gzipped FASTQ files is mostly IO-bound. We could instead focus on making sure we are getting the right answer.

People deride them too often, but this is where a well-written parser is worth it's weight in gold. Heng Li gives us this FASTQ Parser in C.

I downloaded the example tarball and modified the example code (excuse my C...):

#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    int l;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r");
    seq = kseq_init(fp);
    int seqcount = 0;
    long seqlen = 0;
    while ((l = kseq_read(seq)) >= 0) {
        seqcount = seqcount + 1;
        seqlen = seqlen + (long)strlen(seq->seq.s);
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("Number of sequences: %d\n", seqcount);
    printf("Number of bases in sequences: %ld\n", seqlen);
    return 0;
}

Then make and kseq_test foo.fastq.gz.

For my example file (~35m reads of ~75bp) this took:

real    0m49.670s
user    0m49.364s
sys     0m0.304s

Compared with your example:

real    0m43.616s
user    1m35.060s
sys     0m5.240s

Konrad's solution (in my hands):

real    0m39.682s
user    1m11.900s
sys     0m5.112s

(By the way, just zcat-ing the data file to /dev/null):

real    0m38.736s
user    0m38.356s
sys     0m0.308s

So, I get pretty close in speed, but am likely to be more standards compliant. Also this solution gives you more flexibility with what you can do with the data.

And my horrible C can almost certainly be optimised.

Same test, with kseq.h from Github, as suggested in the comments:

My machine is under different load this morning, so I've retested. Wall clock times:

OP: 0m44.813s

Konrad: 0m40.061s

zcat > /dev/null: 0m34.508s

kseq.h (Github): 0m32.909s

So most recent version of kseq.h is faster than simply zcat-ing the file (consistently in my tests...).


What kind of analysis is practically done on GSE data files?

Have you ever used GEO2R? You can use GEO2R to analyze the data dividing it in subgroups (clusters) that interest you.


The effects of incomplete bisulfite conversion upon mapping efficiency

Bisulfite conversion efficiency has no effect on the mapping rate in bismark and similar tools. The reason is that the reads are fully bisulfite converted in silico before alignment to minimize mapping bias.

I would suggest that you play around with the settings handed to bowtie2, such as using local alignment and modifying the --score-min option to allow more mismatches. Alternatively, you might try a different aligner like bwameth, which will always do local alignment.


Why models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decrease

Reading Kaufmann and van Oudenaarden (2007), it seems to validate the first alternative (using results from the Central Limit Theorem):

  Although biochemical fluctuations influence all stages of gene expression, those involving molecules in extremely low abundance are expected from a statistical standpoint to be larger in magnitude and therefore to contribute disproportionately to the overall variation between cells.

Kaufmann, B. B. & van Oudenaarden, A. Stochastic gene expression: from single molecules to the proteome. Curr. Opin. Genet. Dev. 17, 107–112 (2007).


Both the ETE Toolkit and DendroPy offer methods of node removal and reattachment - I would think it would be easy enough to hack together an NNI movement function with them.

You could also investigate the source code of BioPython's NNITreeSearcher() function.


How can I extract normalized read count values from DESeq2 results?

The normalized counts themselves can be accessed with counts(dds, normalized=T).

Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. Given your previous question, we can see that geno_treat has a bunch of levels, which means that expanded models are not in use. In such cases, the baseMean should be the mean of the base factor in geno_treat.


Using the t-SNE algorithm on microarray data + an error bonus

Converting your data.frame to a matrix (and then removing the data.frame) will often free up enough memory that you won't run into this. Note that a matrix is more memory efficient than a data.frame and you're requiring Rtsne() to hold both in memory at the same time (many math-centric functions will end up converting things to a matrix at some point for efficiency).

For what it's worth, it's never been entirely clear to me what the interaction is between a data.frame and the pointer protection stack, but it's often the case that this solves this sort of error.


running a samtools command for multiple bam files from 1000 genomes project

The command you're looking for is:

samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

For a series of intervals, you'd be better off scripting something with pysam:

import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

You might see if there are bigWig files available, they'd be faster to process.


What is a good pipeline for using public domain exomes as controls?

I gather that you want to use the background SNP frequency as a prior for input to your SNP calling algorithm?

I'm not sure of a canned algorithm for doing this, but a quick google shows up some promising links:

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3593722/
https://software.broadinstitute.org/gatk/documentation/article.php?id=4723
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-015-0489-0


Plot information stored in dataframe within a tree (ggtree)

You can accomplish this using the ggtree package available on Bioconductor.
First you will need to combine your tree with the data.
library(tidyverse)
library(ggtree)

ftree <- tree$edge %>%
   as_tibble() %>%
   mutate(Node_number = 1:n()) %>% # finds edge numbering
   right_join(data, by = "Node_number") %>% # find internal node associated with your node numbering
   left_join(fortify(tree, ladderize = FALSE), ., by = c("node" = "V2")) # merge tree and data

Note that this will produce a warning, but it does work.
Then plot the tree using ggtree and its geom_nodepoint layer.
ggtree(ftree, ladderize = FALSE) + # plots tree
   geom_tiplab() + # adds tip labels
   geom_nodepoint(aes(subset = !is.na(Number), size = Number, x = branch)) + # adds your data as points
   scale_size_continuous(range = c(4, 12)) # scales point sizes

You can adjust the size of the points with the scale_size_continuous function.
Data
data <- read_csv('Node_number,Number
50,31
1,15
13,25
35,12
40,4
5,2')

and tree is the same as you provided.
Ggtree also has a decent manual available at https://yulab-smu.top/treedata-book/.


Prediction model for spread of coronavirus

The original R0 calculation was developed in the early 1900s to model malaria transmission and is known as the Ross-MacDonald model so your idea is over 100 years old. The model was very famous for identifying the optimal strategy for malarial control by targeting the adult mosquito rather than the mosquito larvae. What it said was reducing the mosquito population and reducing malaria transmission were two different issues. So the bednet strategy is based in this idea and particularly spraying the inside of houses with insecticide.

R0 was famously used in the 1980s to predict the HIV pandemic when WHO thought it would remain localised in Africa.

The problem with R0 and Coronaviruses is we don't know the asymptomatic population who could also transmit the virus.

This population is potentially very large. Epidemiologists will attempt extrapolation from cruise ship quarrentine, eg the Diamond Princess where everyone is tested to understand this transmission parameter and project this onto the rest of the planet. Clearly this could have limitations or inaccuracies. The other issue I am baffled by is seasonal transmission and whether this will occur with COVID-19.

One thing R0 of COVID-19 will identify is the latent period transmission which occurs before the symptoms, i.e. the incubation period. Essentially, people with cough/mild symptoms before they self isolate and become really ill with a lung infection will transmit the virus. In this context R0 will be/was very influential for understanding when to move between quarrentine and containment phases of the pandemic.


How do you visualize minimap2's PAF output format?

As described in the minimap2 manual, -x STR is used to specify a preset option like -x map-pb for PacBio/Oxford Nanopore read to reference mapping (-Hk19).
In general, I use the -a flag to output alignments in sam format that can then be  sorted and converted to bam format, and viewed in any genome browser. The support for PAF format when it comes to visualization appears to be limited. A few options that I have found, but not tested are:

minidot
pafr
pafpy -- this one may not have the plotting capabilities though


Is it possible to do alignment within Python? Check variants against reference?

As others have said, doing whole-genome (or even whole-chromosome) alignments is the wrong solution. Simply track where you're creating SNVs. If you wrote your SNV locations to a sorted BED file you could use something like the following:

#!/usr/bin/env python
minWidth = 100000  # Only report regions of >=100kb
last = [None, None]
for line in open("some bed file"):
    cols = line.strip().split()
    if cols[0] != last[0]:
        if int(cols[2]) >= minWidth:
            print("{}:1-{}".format(cols[0], cols[2]))
    else:
        if int(cols[2]) - last[1] >= minWidth
            print("{}:{}-{}".format(cols[0], last[1], cols[2]))
    last = [cols[0], int(cols[2])]

This will miss the last part of each chromosome/contig, since it doesn't know how long they are, but you should be able to figure out how to handle that. The memory requirements of this are trivially small.

I suspect there's also a bedtools or bedops command that will directly give you the distance to the next downstream region. One could alternatively parse that (it's all the above code is doing, but using bedtools or bedops would produce more validated results than a quickly thrown together script).


What is the right way of calculating a Phred score by hand?

Not sure if my explanation is any good but let's try...
First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2).
Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550:
n = 10000
(n * 0.01 * 0.5) + (n * 0.1 * 0.5) # -> 550 wrong calls

So the overall error rate is 550/10000 = 0.055 which corresponds to Q=12.6. If you follow your Way 2 you do get this answer:
((n * 0.5 * 0.01) + (n * 0.5 * 0.1))/n = 0.055
-10*log10(0.055) = 12.6

If you follow your Way 1 you get Q=15, which means an overoptimistic error rate of ~0.03 or ~316 wrong calls. Why is that?
Because phred is on a non-linear, logarithmic, scale: a unit change  towards the high values is less important than a unit change towards the lower end. For example, across 10000 bases a change from Q=30 to Q=31 means a change from 10 to ~8 wrong calls - a difference of 2 only. But a change from  Q=10 to Q=11 means 1000 to ~800 wrong calls: the difference is a substantial ~200 calls.
Having said this, the straight average of phred scores is not "wrong". It measures the overall phred score but it does not reflect the expected number of wrong calls, which is probably what you want in this case.


snakemake get environment variable for rule name and log path

You could add an exit trap to the beginning of the shell scripts inside your Snakefile, which is executed after each script is done and sends a message when there was an error:

 shell:
    """
    trap 'if [ $? -neq 0 ]; then ping_slack -m "Error in rule xxx with input {input}"; exit 1' EXIT
    ...


$f in your command is not pointing anywhere since your for loop was defined as for file in *.vcf;. You should use $file instead of $f.

EDIT AS THE OP HAS EDITED THE QUESTION:

First of all, please copy and paste the error messages you are getting when seeking answers from the community, most of the time the error messages are clear enough to point to the right direction.

In you case, the problem is that you are trying to use a variable that you did not declare in the first place. It is explained above in my original answer but here I go again:

With this bit for file in *.vcf;, you are creating a variable file that you will using throughout the iteration. But then within the iteration, the next line (... -operation gx,r,f,f,f "$f" > ..., you try to use another variable file, that has not been declared yet (as far as I can see from your code bit).


Question about public availability of human SNP dataset with country of origin

1000 genomes, HGDP and SGDP all contain SNP data from individuals where the country of origin is known.
They all have slightly different numbers of people from different populations (1000 genomes and HGDP have more individuals in each pop, but fewer pops, and SGDP is the opposite).


pt <- table(Idents(myobj), myobj$orig.ident)
pt <- as.data.frame(pt)
pt$Var1 <- as.character(pt$Var1)

ggplot(pt, aes(x = Var2, y = Freq, fill = Var1)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  scale_fill_manual(values = brewer.pal(12, "Paired")) +
  theme(legend.title = element_blank())

This shall work. myobj is your seurat object. This is my old code so if I remember correct, orig.ident is your treated and untreated conditions. Idents should give you cell clusters.
If this does not work, please let me know


Which taxonomy data type to use for abundance graphs and statistical analyses?

The fundamental answer to this question is the role of allele frequency, i.e. population genetics, within taxonomy. The answer of course is that population genetics and population dynamics should be a integral part of taxonomic designation PROVIDING it is not subject to unreported sample bias. Thus it has to be option 1.

However, your definitions I'm not so sure about because the counts per OTU using your literal definition means every SNP within a species. I don't believe that is what you mean, but that is what you have said. Obviously at a taxonomic focus of a single SNP this isn't  meaningful because there would be as many OTUs as there are sequences and the species definition would be meaningless.

Both approaches suffer from lacking a clear definition to make a species assignment and this is controversial. If you are using a tree to assign species membership, e.g. 16S at 5% identity, then that would generally be acceptable, but some would question it.

You can of course use deep learning algorithms to do this, but this is only as accurate your training and at present (check Bioinformatics journal) the accuracy is at best 90%, which is very low by comparison to phylogeny.

Thus I think it is a non-question, because allelic diversity is extremely important in vast swathes of biological analysis.  Taxonomy is not a disapline in its own right, but a foundational step biologists use to understand the ecotope, ecoregion etc... The only  justification for pursing approach 2 is if there is persistent concern in the bias of the sampling process, to the point where the alleleic (species) frequency becomes meaningless.


How to get the genomic sequences from a blat result?

Via BEDOPS convert2bed (psl2bed) and bedops operations:
$ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed

The file answer.bed will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.
This BED file can be run through samtools faidx or similar to get sequence data.
References:

https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range


`TxDb.Hsapiens.UCSC.hg19.knownGene` and `EnsDb.Hsapiens.v75` can't be installed

I am able to successfully install these packages using:
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))

This completes without any errors / warnings (R v4.0.3 / macOS 10.16), so I don't believe the packages themselves/Bioconductor are 'broken', but I suspect you are having compiler issues due to the Big Sur update.
These are the steps I used to solve the problem:

Reinstall xcode command line tools (even if it says "up to date")

sudo rm -rf /Library/Developer/CommandLineTools
sudo xcode-select --install

Install gcc & llvm via Homebrew (instructions for installing Homebrew) or, if you already have gcc/llvm installed via Homebrew, skip to the next step

# This can take several hours
brew install gcc
brew install llvm

Once you have gcc & llvm installed via Homebrew:

brew cleanup
brew update
brew upgrade
brew reinstall gcc
brew reinstall llvm

Link some headers into /usr/local/include

sudo ln -s /Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include/* /usr/local/include/

# You can safely ignore warnings like this:
#ln: /usr/local/include//tcl.h: File exists
#ln: /usr/local/include//tclDecls.h: File exists
#ln: /usr/local/include//tclPlatDecls.h: File exists
#ln: /usr/local/include//tclTomMath.h: File exists
#ln: /usr/local/include//tclTomMathDecls.h: File exists
#ln: /usr/local/include//tk.h: File exists
#ln: /usr/local/include//tkDecls.h: File exists
#ln: /usr/local/include//tkPlatDecls.h: File exists

Create a new ~/.R/Makevars file (if you already have a ~/.R/Makevars file, save it in a different directory (away from ~/.R/)) and include only these lines in the file:

FLIBS=-L/usr/local/gfortran/lib/gcc/x86_64-apple-darwin19/10.2.0 -L/usr/local/gfortran/lib -lgfortran -lquadmath -lm
CXX1X=/usr/local/gfortran/bin/g++
CXX98=/usr/local/gfortran/bin/g++
CXX11=/usr/local/gfortran/bin/g++
CXX14=/usr/local/gfortran/bin/g++
CXX17=/usr/local/gfortran/bin/g++

LLVM_LOC = /usr/local/opt/llvm
CC=/usr/local/gfortran/bin/gcc -fopenmp
CXX=/usr/local/gfortran/bin/g++ -fopenmp
CFLAGS=-g -O3 -Wall -pedantic -std=gnu99 -mtune=native -pipe
CXXFLAGS=-g -O3 -Wall -pedantic -std=c++11 -mtune=native -pipe
LDFLAGS=-L/usr/local/opt/gettext/lib -L$(LLVM_LOC)/lib -Wl,-rpath,$(LLVM_LOC)/lib
CPPFLAGS=-I/usr/local/opt/gettext/include -I$(LLVM_LOC)/include -I/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk/usr/include

Compile a package from source in R/Rstudio

# To check whether openmp is enabled, compile data.table:
install.packages("data.table", type = "source")

If your package fails to compile, a couple of SO users have had to install a fresh gfortran (re: https://stackoverflow.com/a/65334247/12957340), which you can download from https://github.com/fxcoudert/gfortran-for-macOS/releases/tag/10.2-bigsur-intel

Once you can install packages from source without error/issue, update Bioconductor
BiocManager::install(version = "3.12")

Then, finally, you should be able to install your packages using
BiocManager::install(c("Homo.sapiens", "EnsDb.Hsapiens.v75", "TxDb.Hsapiens.UCSC.hg19.knownGene"))


problem of "ordering in physical positions" phasing SNPs with Shapeit

The manual suggests that each "chromosome" needs to have its own input file.

You have to split the dataset by chromosomes. ShapeIT can just phase one chromosome at a time. If you are working on GWA data, you can proceed as follows:

for chr in $(seq 1 22) ; do plink --file example --chr $chr --recode --out example_chr$chr ; done

You will obtain 44 files; example_chr1.ped, example_chr1.map, ..., example_chr22.ped, example_chr22.map that you can phase separately.


Finding out proteins that have a similar domain to a specific toxin

A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.

This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast results do include conserved domain information though.


Normalization methods with RNA-Seq ERCC spike in?

You may consider using RUVSeq. Here is an excerpt from the 2013 Nature Biotechnology publication:

  We evaluate the performance of the External RNA Control Consortium (ERCC) spike-in controls and investigate the possibility of using them directly for normalization. We show that the spike-ins are not reliable enough to be used in standard global-scaling or regression-based normalization procedures. We propose a normalization strategy, called remove unwanted variation (RUV), that adjusts for nuisance technical effects by performing factor analysis on suitable sets of control genes (e.g., ERCC spike-ins) or samples (e.g., replicate libraries).

RUVSeq essentially fits a generalized linear model (GLM) to the expression data, where your expression matrix $Y$ is a $m$ by $n$ matrix, where $m$ is the number of samples and $n$ the number of genes. The model boils down to

$Y = X*\beta + Z*\gamma + W*\alpha + \epsilon$

where $X$ describes the conditions of interest (e.g., treatment vs. control), $Z$ describes observed covariates (e.g., gender) and $W$ describes unobserved covariates (e.g., batch, temperature, lab). $\beta$, $\gamma$ and $\alpha$ are parameter matrices which record the contribution of $X$, $Z$ and $W$, and $\epsilon$ is random noise. For subset of carefully selected genes (e.g., ERCC spike-ins, housekeeping genes, or technical replicates) we can assume that $X$ and $Z$ are zero, and find $W$ - the "unwanted variation" in your sample.


Understanding the initialization of the Needleman-Wunsch algorithm

The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment

The dynamic programming is useful because you need first to calculate the score for each position, and then once you have the score, you can calculate the optimum path. So you solve your problem by first calculating what would happen if (there was a gap, there was a mismatch...) and then you find the optimum alignment.


How to estimate whether a long-read is meaningful sequence?

As a first pass, you could check if the read is chimeric. Porechop searches for known nanopore adaptors both on the ends and through the middle of the read. This won't resolve issues around blocked or empty pores, but it will at least check if you have found two long-ish reads lumped into the same file.

By default, Porechop splits chimeric reads into two (or more, I suppose!) non-chimeric reads, but the option --discard-middle would be a quick and easy way to check - run it on a fasta file containing only your long read, and if the output is empty, the read is chimeric.


RNASeq analysis using featureCount and EdgeR

It makes no difference if you process the BAM files one at a time with featureCounts or all together, except that it changes how you have to read the files into R.
You can supply edgeR with lists of contrasts to have it compute fold-changes and p-values for. Please have a look at the edgeR user guide for examples.


TagReadWithGene missing when using latest version of Drop-seq_tools

As per this response from Drop-seq_tools team, at the moment intronic regions cannot be annotated and TagReadWithGeneExon should be used instead.

However, in the version of the public repository:

  TagReadWithGeneExonFunction is the old code, for people who want the more simple tagging methodology and rely on it for 3rd party tools, and TagReadWithGeneFunction, which for your purposes replaces TagReadWithGeneExon. [...]  You’ll see additional options to extract expression for CODING, UTR, INTRONIC reads [...]


Estimating BAM file from compressed fastq file size

A BAM file should have a roughly similar size to a compressed FASTQ file, because they're both compressed files that contain similar information. BAM will be slightly larger, because it contains information about mapping and differences from the reference.

If you are concerned about space and have a reference sequence available, you can store mapped reads as CRAM format, which will save a lot more space.


What is local realignment and what is the problem it solves?

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get:

ACTGACTGACTGTac    (lower case denotes soft-clipping)
ACTGACTGACTGTTAAGAACGACTACGAC

or

ACTGACTGACTGTAC    (N.B., you've created some possible SNPS)
ACTGACTGACTGTTAAGAACGACTACGAC

But what happens if you look at a bunch of reads? See below:

ACTGACTGACTGTAC
  TGACTGACTGTACTA
    ACTGACTGTACTAAG
ACTGACTGACTGT--TAAGAACGACTACGAC

I've made the 2 base insertion obvious. As you can see, using many alignments, we can "realign" our original read to see that it actually ends in a 2 base insertion, something that's completely impossible to find with a single read. This realignment is part of what the GATK haplotype caller does (it makes a local de Bruijn graph). That's also what my MethIndelRealigner does. I have another real-world example in its documentation.

Raw alignments in IGV:

Alignments after realignment:

One can find plenty of similar scenarios that result in SNPs rather than simply losing coverage for InDels.


Low percentage of reads with consistent barcodes in Split-Seq run

In my case the problem was due a wrong set of allowed round1 barcodes, therefore the filtering pipeline was removing most part of the reads. In particular, with the right barcode1 whitelist the total number of reads after filtering is 17M (about 100 times more).


It depends on what you're looking for. If you're just looking for sequence homology, then you can simply pick the best hits from a blast search. If, however, you are referring to functional homology, if you are looking for the protein which has the same functions as your query, then it's more complicated.

Sequence homology is not enough to infer functional homology. For example, you can have cases of gene duplication and subsequent functional divergence. Such paralogs are still homologs (paralogs are a subset of homologs), but they don't necessarily have the same function. It is also often the case that the homolog (be it orthologous or paralogous) of a protein in species B has a completely different function than its homolog in species A despite a high level of sequence similarity. This is usually very hard to determine in silico.

To find the functionally homologous protein(s), you would ideally need to identify the essential residues that allow your protein to perform its function. This could be done using something like PFam which will identify protein domains. You can then check whether the homologs you find also have this domain.

This is essentially what PSI-blast does. Although it doesn't take domains into account, each successive iteration is used to build a model of your proptein. The model is built under the reasonable assumption that highly conserved residues are important. So it will consider more diverged sequences as homologous if those residues are conserved.

If you know how your protein works and what residues are important, you can use that knowledge to refine the results of your PSI-blast. If you don't, you'll have to use only "good" hits to make the model. One way to do this, for well studied proteins, is to only add proteins that are already annotated as homologs of what you are looking for to build your model, then use that model to search in un-annotated species.

If you don't know, you could try looking for recognizable protein domains in your query protein (use PFam) and then use the HMM (hidden markov model) of the domain to identify important residues. For example, this is the HMM logo for the zf-A20 zinc finger domain:

The huge cysteine (C) residues are shown that size because they are very conserved across proteins carrying this domain and, presumably, are functionally important for the domain. So, if you pass your protein through PFam and identify domains, find the important residues and make sure all your homologs have those conserved. If using PSI-blast, only include sequences where those residues are conserved in the results you keep.

Finally, another useful tool that works in the same way is HMMER. This takes a protein alignment as input, like PSI-blast builds an HMM model from it and then can use this model to query a protein database for more hits. Methods like HMMER and PSI-blast are far better than simple sequence similarity approaches when looking for homologs.


How to use the DI-tector software for defective genomes analysis?

See the help:

python DI-tector_06.py -h
usage: DI-tector_06.py [-h] [-g HOST_REF] [-s MIN_SEGMENT] [-m MIN_MAPQ]
                       [-n NB_READS] [-o OUTPUT_DIRECTORY] [-t TAG] [-d]
                       [-l INDEL_LENGTH] [-f] [-p POLARITY] [-q] [-k]
                       [-x NB_THREADS]
                       Virus_Ref Input_Data

positional arguments:
  Virus_Ref             Virus genome reference sequence in FASTA format.
  Input_Data            File containing single reads in FASTQ format.

optional arguments:
  -h, --help            show this help message and exit
  -g HOST_REF, --Host_Ref HOST_REF
                        Host genome reference sequence in FASTA format.
  -s MIN_SEGMENT, --Min_Segment MIN_SEGMENT
                        Minimum segment length. Default is 15.
  -m MIN_MAPQ, --Min_MAPQ MIN_MAPQ
                        Skip alignments with MAPQ smaller than INT. Default is
                        25.
  -n NB_READS, --Nb_Reads NB_READS
                        Show only DVGs with counts reads > or egal to INT.
                        Default is 2.))
  -o OUTPUT_DIRECTORY, --Output_Directory OUTPUT_DIRECTORY
                        Enter a directory name that all compiled output files
                        will be saved in.
  -t TAG, --Tag TAG     Enter a tag name that will be appended before each
                        output file. Default is 'DI-tector'.
  -d, --DVG_sequences   Generate multi-fasta file with DVG sequences. Default
                        is (OFF).
  -l INDEL_LENGTH, --InDel_Length INDEL_LENGTH
                        Skip alignments with size of InDels smaller or egal to
                        INT. Default size is 1.
  -f, --Fasta           Select '-f' if data is in FASTA format fasta. Default
                        is FASTQ.
  -p POLARITY, --Polarity POLARITY
                        [0] Positive strand genome / [1] Negative strand
                        genome. Default is 0.
  -q, --No_Quantification
                        Inactive percentage quantification. Quantification
                        need bedtools Default is (ON).
  -k, --Keep_files      Keep intermediaite files (i.e. alignment etc...).
                        Default is (OFF).
  -x NB_THREADS, --Nb_threads NB_THREADS
                        Number of threads. Default is 1.


Your table is not very clear. If I understood it correctly, there are three groups of t-cells, SPtetramer+, Ly49- and Ly49+ and within it there is a MOGSP vs MOG comparison. I would set colData like this:
                   name             Ly treatment
1  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
2  SPtetramer+CD8+TCell SPtetramerplus     MOGSP
3  SPtetramer+CD8+TCell SPtetramerplus       MOG
4  SPtetramer+CD8+TCell SPtetramerplus       MOG
5  SPtetramer+CD8+TCell SPtetramerplus       MOG
6        Ly49-CD8+TCell      Ly49minus     MOGSP
7        Ly49+CD8+TCell       Ly49plus     MOGSP
8        Ly49-CD8+TCell      Ly49minus     MOGSP
9        Ly49+CD8+TCell       Ly49plus     MOGSP
10       Ly49-CD8+TCell      Ly49minus     MOGSP
11       Ly49+CD8+TCell       Ly49plus     MOGSP

I generate an example dataset:
counts = matrix(rnbinom(11000,mu=100,size=1),ncol=11)

And in the design you specify these two factors, meaning gene expression can be affected by these two conditions:
dds = DESeqDataSetFromMatrix(counts,colData=colData,~Ly +treatment)
dds = DESeq(dds)

You can compare them like this:
head(results(dds,c("Ly","Ly49plus","Ly49minus")))
log2 fold change (MLE): Ly Ly49plus vs Ly49minus
Wald test p-value: Ly Ly49plus vs Ly49minus
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915  0.712987165797688 1.28291188260395  0.555756927241586
2 81.7945234843954 -0.418266530066077 1.14378290088875 -0.365686993345568
3 94.4351681904941 -0.775592496878288 1.32965355129901 -0.583304196886903
4 107.484335203298  -1.36694885973712 1.38595911743023 -0.986283680770933
5 61.2205696982841 -0.760219864564472 1.16779888688769 -0.650985262188887
6 95.5908168437137  0.497427543096335 1.12768919469662   0.44110340458672
             pvalue              padj
          <numeric>         <numeric>
1 0.578377034355333 0.977537706056193
2    0.714598652706 0.992321035598301
3 0.559688537757548 0.977537706056193
4 0.323993925580938 0.954772990635094
5 0.515056000552832 0.977537706056193
6 0.659138138867294 0.992321035598301

To compare groups within another factor:
results(dds,c("treatment","MOGSP","MOG"))
log2 fold change (MLE): treatment MOGSP vs MOG
Wald test p-value: treatment MOGSP vs MOG
DataFrame with 6 rows and 6 columns
          baseMean     log2FoldChange            lfcSE               stat
         <numeric>          <numeric>        <numeric>          <numeric>
1 91.4458114779915 -0.875905107403943 1.43007567177694 -0.612488642867122
2 81.7945234843954  0.595588650280833 1.28109080878515  0.464907441530728
3 94.4351681904941 0.0255578400583986 1.48592401497315 0.0171999643325371
4 107.484335203298  -1.46652057431159 1.55285829334496 -0.944400774105796
5 61.2205696982841  -1.27442596494316 1.30489822606005  -0.97664778715433
6 95.5908168437137  -2.35402796630456 1.26608624796739  -1.85929510733078
              pvalue              padj
           <numeric>         <numeric>
1  0.540214509839105 0.944881849859735
2  0.641997741235223 0.973651054534338
3  0.986277090644156 0.994944323499595
4  0.344964886056295  0.87421063696937
5  0.328743552427222                NA
6 0.0629853201285643 0.525889480109579

The colData i used:
structure(list(name = structure(c(3L, 3L, 3L, 3L, 3L, 1L, 2L,
1L, 2L, 1L, 2L), .Label = c("Ly49-CD8+TCell", "Ly49+CD8+TCell",
"SPtetramer+CD8+TCell"), class = "factor"), Ly = c("SPtetramerplus",
"SPtetramerplus", "SPtetramerplus", "SPtetramerplus", "SPtetramerplus",
"Ly49minus", "Ly49plus", "Ly49minus", "Ly49plus", "Ly49minus",
"Ly49plus"), treatment = structure(c(2L, 2L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("MOG", "MOGSP"), class = "factor")), class = "data.frame", row.names = c(NA,
-11L))


How to edit the headers of multiple fasta files from multiple folders?

You were so close! This should do the trick:

for i in {314..727}; do
    sed -i.bak "s/^>/>${i}_/" mgm4761${i}.3/*.fna
done

Note that I am using -i which makes sed edit the original file. Because I use -i.bak, this will create a backup file with the original data. Once you are sure the command does what you need, you can delete all of these files with rm */*fna.bak.


What is the `table` parameter in the UCSU ldHgGene tool?

I stumbled across this issue after getting the same warning message from liftOver. By using the -out option of lgHgGene, the database parameter is ignored but must still be provided. The table parameter should be the second column from your GTF file. Note that you must also specify the -gtf option (otherwise it expects GFF format). In my case, the GTF file looks like:

chr1    mm10_knownGene  exon    3073253 3074322 0.000000        +       .       gene_id "ENSMUST00000193812.1"; transcript_id "ENSMUST00000193812.1";
chr1    mm10_knownGene  exon    3102016 3102125 0.000000        +       .       gene_id "ENSMUST00000082908.1"; transcript_id "ENSMUST00000082908.1";
chr1    mm10_knownGene  exon    3206523 3207317 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3213439 3215632 0.000000        -       .       gene_id "ENSMUST00000159265.1"; transcript_id "ENSMUST00000159265.1";
chr1    mm10_knownGene  exon    3205901 3207317 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
chr1    mm10_knownGene  exon    3213609 3216344 0.000000        -       .       gene_id "ENSMUST00000162897.1"; transcript_id "ENSMUST00000162897.1";
...

By default, ldHgGene outputs a leading 'bin' column which won't work with liftOver; this can be avoided with the -nobin option. So my command was:

ldHgGene -gtf -nobin -out=mm10.GENCODE_VM20.gp ignored mm10_knownGene mm10.GENCODE_VM20.gtf

The resulting genePred output can be fed directly into liftOver -genePred:

ENSMUST00000193812.1    chr1    +       3073252 3074322 3074322 3074322 1       3073252,        3074322,
ENSMUST00000082908.1    chr1    +       3102015 3102125 3102125 3102125 1       3102015,        3102125,
ENSMUST00000162897.1    chr1    -       3205900 3216344 3216344 3216344 2       3205900,3213608,        3207317,3216344,
ENSMUST00000159265.1    chr1    -       3206522 3215632 3215632 3215632 2       3206522,3213438,        3207317,3215632,


Depending on the cluster management tool, you might have received e-mails when the "job" begins and ends. If so, you can check the "Exit status" of the job. For example, in the case of our HPC the relevant lines from the "job completion e-mail" would be:

Execution terminated
Exit_status=0

BTW, at first I was confused with the phrase "Execution terminated", however, in the scheduler used in our system, it means "the run has been successfully or unsuccessfully completed". Exit_status=0 int the next line means a successful run and for the unsuccessful runs there is a list with the error codes.

Secondly, if your run was successful, you would see new fastq files containing "_trimmed" in the filename (reading this from your code). These fastqs should be more or less similar in size to the original depending on how aggressive the trimming settings are.

If you can provide more information, i.e. the type of cluster management tool used, the log file, ..., you can get more help and precise answers from the community.


Finding SNP/Indel Markers of Interspecific Variation

Problem statement
You want to find a set of SNPs in populations of two species, for each of which you have a distinct reference genome. You then want to know which of those SNPs are shared between the two genomes (there is variation in each species), and which ones are private to only one of the species (in the other genome there is no variant at that location).
To accomplish this task, you need to first find SNP variation in each genome, and then figure out which SNPs are the same between the two genomes. In other words, line up two VCF files for different genomes and figure out which SNPs are the same between them. In comparing genome assemblies, this is sometimes called a liftover operation.
Genomes can have different coordinate systems.
Let's first think about what exactly you are looking at with 2 different reference genomes. For example, here are outputs of samtools faidx for 2 different S. cerevisiae strains, S288C and Kyokai-7:
# S288C
ref|NC_001133|  230203  95  60  61
CHR2.19980913   813142  234173  60  61
CHR3.19980521   315339  1060907 60  61
CHR4.19990210   1531929 1381540 60  61
CHR5.19970727   576870  2939039 60  61
CHR6.19970727   270148  3525562 60  61
CHR7.19970703   1090936 3800252 60  61
CHR8.19970727   562638  4909411 60  61
ref|NC_001141|  439885  5481523 60  61
CHR10.19970727  745440  5928778 60  61
CHR11.19970727  666448  6686681 60  61
CHR12.19970730  1078172 7364277 60  61
ref|NC_001145|  924430  8460517 60  61
ref|NC_001146|  784328  9400452 60  61
CHR15.19970811  1091283 10197892    60  61
ref|NC_001148|  948061  11307461    60  61
ref|NC_001224|  85779   12271444    60  61

# Kyokai first 10 lines
gi|347729946|dbj|BABQ01000001.1|    5491    137 80  81
gi|347729945|dbj|BABQ01000002.1|    4704    5834    80  81
gi|347729944|dbj|BABQ01000003.1|    203046  10734   80  81
gi|347729943|dbj|BABQ01000004.1|    1036    216456  80  81
gi|347729942|dbj|BABQ01000005.1|    126714  217642  80  81
gi|347729941|dbj|BABQ01000006.1|    121900  346077  80  81
gi|347729940|dbj|BABQ01000007.1|    20385   469638  80  81
gi|347729939|dbj|BABQ01000008.1|    243458  490415  80  81
gi|347729938|dbj|BABQ01000009.1|    189876  737054  80  81
gi|347729937|dbj|BABQ01000010.1|    81973   929441  80  81
# ...
# goes on for many thousands of lines

What you will notice is that there are significant differences between the 2 assemblies. The first is almost all in chromosomes, the second is largely mapped to various chromosomes but exists in thousands of tiny contigs. There is no way to just map reads to them, call variants, and assume that you (or bcftools) can figure out what is going on.
This is possibly an extreme example, but even in the case of two genomes assembled to chromosome pseudomolecules, there is basically no chance that those two chromosomes are on the same coordinate system or have similar metadata. That means again that bcftools probably won't be able to figure it out, because a SNP at position 112929 in genome A will be homologous to position 110201 in genome B.
Liftover: a way to solve the mapping problem between 2 genomes.
You can create a mapping between the two genomes. This is 100% a solved problem in bioinformatics, and it is more or less the same as doing liftovers between different assembly versions of the same genome. I think that you can more or less just use the UCSC LiftOver tool to convert SNP coordinates, unless I am missing something. That how-to guide that I linked is going to give you much better instructions than I can give you, here is a more minimal version that might be simpler.
The only extra suggestion that I have is that the guide uses BLAT to align the two assemblies, and minimap2 is a much faster tool. You may need to fight the input formatting a little bit, but if you run into runtime issues you might make that change. This is why I suggested it in the comment, though I wasn't aware that the liftover docs were so complete.
Short pipeline:

Align reads (individuals of species A to reference A, individuals of species B to reference B).
call variants-->get VCF files using preferred workflow for each of A and B independently.
compute liftover mapping (chain file, usually) from reference A to reference B.
use computed liftover mapping to actually lift VCF over from one genome to another.
align A and B references, call SNPs between them (possibly can be done from (3) output). Ensure that you are using as reference the liftover target genome so that the VCF is informative.
For each SNP private to each species (output of 5), query the BAMs to ensure that no variation exists.


Using Python, how to convert a pandas DataFrame into a VCF?

You will have to define the header either from an existing VCF file or hardcoded into you python script.  Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.

You will need to order the dataframe columns to match the VCF spec and insert the empty/ default ID, FILTER, INFO and QUAL columns

import pandas as pd

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']]

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

Then myfile.vcf looks like

##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40
20  17330       T   A   40
20  1110696     A   G   40
20  1230237     T   C   40
20  1234567     GTCT    G   40


Which software is used to run Molecular Dynamics simulation of a small RNA hairpin?

Google found this review that might be of interest: https://pubs.acs.org/doi/10.1021/acs.chemrev.7b00427

Here is a link to a free full-text version:
https://www.ncbi.nlm.nih.gov/pmc/articles/pmid/29297679/

The authors "here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA–ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field."

Hope this helps.


Loss of predictive power of polygenic risk score when dataset contains missing variants

Which software package around you using, PRSice? Generally missing values comprising half the total observations are extremely difficult to calculate using regression analysis (or many statistical analyses for that matter.
One way forward is to include (an) example(s) of all 40 variants into your data set and perform the analysis with and without the added data, then compare this against a data set where all 40 variants are observed and your samples are excluded. You might find the missing values have low probability of PRS and this would then explain your results.
Missing analysis in regression is a specialist area of statistics and is not for the faint hearted, so a work-around is the best way forward.


Interpreting GWAS results with different settings

Based on the QQ plots, filter out variants with p > 0.01, then display the beta statistic (or whatever other statistic you're using to demonstrate association). If you want to incorporate the p-value in the displayed results, show the minimum test statistic from a 95% confidence interval. Here's an attempt I made at doing that with some GWAS summary statistics relating to household income (something that wouldn't be expected to have a heritable component):
https://twitter.com/gringene_bio/status/1207617723586371584
After looking at the test statistic, you may be able to make better decisions about thresholding with the other values. One thing you could try is to incorporate those values into a scatter plot that shows more than two dimensions. The other statistics could be represented as size, colour, or Z location on a 3D scatter plot.
P-values should not be overemphasised, used for ranking, or displayed on their own. See the American Statistical Association's report on the p-value for more information:
https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i30
As far as I'm aware, the p-values as calculated by GWAS software represent uncertainty in the test statistic (i.e. t, or beta), and not the probability of association. If you want to calculate that association probability, carry out bootstrap sub-sampling of the individuals and then determine the proportion of sub-samples in which the tested variant appears in, say, the highest 5% of variants, ranked by the test statistic. Instead of filtering on HWE, MAF and MCR, you could try filtering on the bootstrap proportion instead, which should incorporate all statistics that introduce uncertainty into the dataset.
More details and figures can be found on this poster:
https://f1000research.com/posters/5-2190
If you're going to be generating a PCA on these data, please consider labelling the cases and controls (or low/high for quantitative traits) on the PCA to verify that the groups look similar in terms of overall population structure. Hopefully you won't find any blobs in the PCA that have only one of the two groups in them, but if you do, they should be dealt with first. Any differences in population structure (even very small, subtle differences) will dominate the GWAS results because structure-associated signal will cover a huge region of the genome in comparison to effects from a few (or a few hundred) SNPs. Here are some papers I've found insightful regarding this issue:

https://doi.org/10.1006/tpbi.2001.1543
https://doi.org/10.1038/s41467-018-05257-7


Many thanks for raising this issue. The web service is now working. Only the glyph related web service calls were affected due to restructuring folders at the SBOL-Visual GitHub repository. I now revived the glyphs using the sbol-visual-ontology GitHub repository. Over the summer, we will develop the web service further so that you can request different versions of glyphs. Currently, the web service works with the SBOL Visual v2.3 glyphs and we will provide support for SBOL Visual v3 glyphs soon.
The web service is available at:
https://vows.sbolstandard.org/
Please let me know if there is anything I can clarify.


changing color key range to specified range in heatmap.2 function

You just need to add more colors. The following, for example, will work:

heatmap.2(d_matrix, scale="none", trace="none", Rowv=F, Colv=F, breaks=seq(0,1,0.01), dendrogram="none", col=redblue(100))

Note that length(seq(0, 1, 0.01)) is 101, therefore redblue() must return one fewer (100) colors (or colours if you prefer extra vowels).


How to find Genotype, from list of mutations:lineage list?

import re
import os
import pandas as pd
import logging
import sys
import csv
import itertools

mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

        combinations = ['_'.join(sorted(x)) for x in combinations]
        combinations = [x.split('_') for x in list(set(combinations))]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

    df = pd.read_csv('Results/' + lineage).dropna()

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

            out_df = pd.DataFrame({'combination':combination_labels,
                               'count':combination_counts})
            out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
            out_df = out_df.sort_values('percentage', ascending = False)

            out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
                     header = True,
                     index = False)


Looking at the htslib source code, it seems like the answer is no, a single region specifier cannot span multliple chromosomes.

Seems like the closest you can get is to enumerate all the chromosomes as separate regions on the command line. So to get everything between chrA:X and chrN:Y, you would have to write:

tabix inputfile.gz chrA:X chrB chrC ... chrN:1-Y

(If you omit the end coordinate like chrA:X, the region encompasses chrA from the specified coordinate to the end. If only specify the chromosome name without the coordinates, like chrB, the region encompasses the whole named chromosome, and for the last chromosome, you can specify up to the desired point by starting from the beginning: chrN:1-Y.)


Cufflinks Error: sort order of reads in BAMs must be the same

Just add --dta-cufflinks to your Hisat2 command, so that the output alignment file provides the attributes needed by Cufflinks (XS flags). This should do the job.

From manual:

  Report alignments tailored specifically for Cufflinks. In addition to
  what HISAT2 does with the above option (--dta), With this option,
  HISAT2 looks for novel splice sites with three signals (GT/AG, GC/AG,
  AT/AC), but all user-provided splice sites are used irrespective of
  their signals. HISAT2 produces an optional field, XS:A:[+-], for every
  spliced alignment.


How do I Investigate and test the hypothesis on the effect of street-light regime on insects?

The hypothesis can be examined by t-test. You can group the data (number of catches) into  dimmed light, and fully lit street lights and proceed with t-test to assess whether dimmed light results  in less en-catchment or higher en-catchment or equivalent encatchment.  Alternatively, a linear model can be applied subject to certain conditions  here. PCA is not  appropriate here.


How is PCR duplication rate computed in scATAC-seq?

The answer is in the supplementary methods page 3.

  Finally, all fragments in the same library with duplicate start and
  end coordinates were removed using Picard


How to import a phylogenetic tree (in nexus format) with node support values in R

I recommend you this tutorial about ggtree. You will find the solution.
[https://guangchuangyu.github.io/ggtree-book/short-introduction-to-r.html][ggtree]


When performing differential expression analysis, should genes with low read counts be removed before or after normalization?

The more genes you have the more robust the scaling factor is (among the reason why one doesn't normalize to ERCC spike-ins without a compelling reason), so I suppose in theory it's better to filter after determining the scale factor. Having said that, I'd be surprised if the results changed much either way. Unless you end up filtering out a LOT of genes the scaling factors should be fairly robust to the change.


In-sample and across samples normalized expression

Yes, this is a standard way of obtaining RPKM/FPKM/CPM values for plotting. Not that you do not need to use a for loop for any of the computations in R. You have a matrix of normalized values and things like transcript length are constant across samples (at least unless you're using something like salmon...although then you'd have TPMs to begin with).


File format of substitution matrix in clustalw

This example seems to have information on this, looks like it should have this format:
   A   G   C   T   U   *
A  1   -1  -1  -1  -1  -1
G  -1  1   -1  -1  -1  -1
C  -1  -1  1   -1  -1  -1
T  -1  -1  -1  1   -1  -1
U  -1  -1  -1  -1  1   -1
*  -1  -1  -1  -1  -1  -1

However, the page I linked is specifically about a bug regarding ClustalW ignoring custom made matrices. It's a pretty old post so it might have been fixed since. However, if you encounter that bug, that same page thankfully does contain a workaround!


The ‹msa› package on Bioconductor does exactly that.

It doesn’t hand the results on a silver platter, though: you’ll need to read the vignette carefully to learn how to use it. But after that it’s pretty powerful.


Counting number of possible alignments between two DNA sequences using python

A closed-form solution is offered in An exact formula for the number of alignments
between two DNA sequences by Torres, Cabada, and Nieto:

$$f(m,n)=\sum_{k=0}^{min(m,n)}2^{k}\binom{m}{k}\binom{n}{k}$$

If this solution seems reasonable, you could calculate this without BioPython, but with a simple xrange loop, the power operator, and scipy.special.binom:

#!/usr/bin/env python

import sys
import scipy.special

def alignments(m, n):
    if m < 0 or n < 0:
        raise ValueError('m and n should be non-negative')
    s = 0
    b = min(m, n) + 1
    for k in xrange(0, b):
        s += (2 ** k) * scipy.special.binom(m, k) * scipy.special.binom(n, k)
    return s

def main():
    m = 4
    n = 2
    a = alignments(m, n)
    assert (m == 4 and n == 2 and a == 41), "something went wrong!"
    sys.stdout.write("alignments(%d, %d) = %d\n" % (m, n, a))

if __name__ == "__main__":
    main()

I don't know what your m and n are, but keep in mind that this solution grows quickly. Also, I did some digging and scipy approximates the calculation of the binomial coefficient, which is probably going to be an issue for large m and n.

If you're working with large values, you may want to look into using the Python long type and consider the accuracy of the answer. The Stack Overflow link offers alternatives to scipy that generate an exact calculation of the binomial coefficient, at the cost of a longer runtime.


How to determine a protein's cellular location based on its sequence?

You could use Blast2go to find the predicted relevant gene ontologies of your sequences. If you look for the cellular compartment sub ontology you should find the location of such protein. I have never used this tool, but here is the paper describing how it works.

I wouldn't restrict to the first predicted match of blastp see if some more genes with a good fit are also located at the same sites. But otherwise your workflow seems right to me.

But to prove this is in X or Y you, or your team or collaborators, should always check with histochemistry if it really allocates on the predicted place. All the softwares predictions must be validated with some biological experiments to prove them.


What are some of the major questions pursued by the field bioinformatics?

Welcome! Bioinformatics is just a term that describes an approach to doing biology. People who use ( DNA sequencing/genomics/protein informatics + programming ) to answer biological questions are after the same answers as the the people who use ( pipettes/microscopes/molecular biology + a lab notebook ).

My best suggestion is to learn about some biological concepts that are interesting to you, see how people have applied bioinformatics to that branch of biology, and think if there is something you might be interested in doing to work on that question. It is critical to have a strong understanding of the underlying biology to get to the core of new findings using bioinformatics.

A lot of the hardcore algorithm development going on today is in genome assembly with long + noisy DNA sequencing reads, aligning long+noisy reads to genomes and trying to figure out the differences between the maternal and paternal haplotypes, working on evolutionary questions with population-level genetic variation data, and figuring out how things are related to each other with phylogenetics, to name a few.

As far as a machine-learning type approach which you seem to be hinting at in your question, one very hot topic at the moment is single-cell sequencing. This is a set of technologies that allows us to get DNA, RNA, or protein sequencing from a population of single cells to look at heterogeneity in tissues, timepoints, et cetera.

Good luck! If you have some free time on your hand this website, ROSALIND, is a fun project Euler-type website where you can learn some biology while programming bioinformatics solutions.


How do I rewrite a read group using pysam?

If the goal is to assign all alignments to a single read group (which this code seems to do), then Picard AddOrReplaceReadGroups might help: https://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups.

Modifying headers using PySAM might be tricky – the internal representation of headers has recently been changed and seems not to be stable across PySAM versions.


Nanopore data clarification regarding merging of samples

It is entirely unclear which type of analysis you aim to do. You probably should be more specific when asking questions.
That said, given that the run was going perfectly fine until a power issue stopped it, all data produced should be fine and you can indeed combine those runs.
You probably shouldn't bother with poretools. If you are using fastq files to start your analysis you can just cat the basecalled files together. Also after alignment, you can merge the bam files. If your analysis requires the fast5 files then you can probably just put these in the same folder, but you should tell us which analysis you are doing.


Fastq: how can I check if they are from DNA or RNAseq data?

There are a few quick'n'dirty ways depending on the type of data.
In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC
which calcualtes the fraction of reads aligning to exon, intron, intergenic etc.
RNA-seq, if you use a standard aligner such as bowtie2 or bwa, and you have a higher species that splices its RNA, then you should see most reads aligning to exons, and quite a fraction of reads that span a splice junction being unmapped (because the exon-intron-exon gap is large and non splice-aware tools do not "bridge" this without tremendous drop in mapping quality so the reads goes unmapped).
DNA-seq (if whole genome) should have a somewhat even coverage across the genome on the global scale. Targeted assays such as ATAC-seq should give distinct peaks that you can easily see in a genome browser just by eye.
If it is exome-seq (so targeted DNA-seq via exon pulldown) then it might get tricky because also most reads (like in RNA-seq) should align to exons, but you should not have spliced reads so the 5' part of the read starting in exon1 and 3' part ending in exon2. Instead many reads should also cover the exon-intron boundary.
As said, I would simply align it and check on a browser, I could imagine one can make a sophisticated guess by just looking at the alignments by eye.
A different topic is reliability, because I personally would not touch a dataset with unknown origin, because that also means you have no idea about the protocols that were used for library prep and what the sample actually is, so I wonder which question your analysis is going to answer.


You're assuming 2-line FASTA format, but the sequence can (and typically does) span multiple lines. You can use Biopython to build a simple parser:

from Bio import SeqIO

for record in SeqIO.parse("file_name.fa", "fasta"):
    if "16S" in record.description:
        print(record.format("fasta"))


Genomic coordinates of predicted binding pocket

It turns out there is a tool that does this beautifully.  I used the proteinToGenome function from the ensembldb package.

The code looks like this:

# Define the Ensemble database object for the
edbX <- filter(EnsDb.Hsapiens.v75, filter = ~ seq_name == "X")
# Set the first and last amino acid you want to convert and the ensembl protein id
GAGE10_prt <- IRanges(start = 5, end = 9, names = "ENSP00000385415")
# Perform the conversion.  GAGE10_prt is an object that has the coordinates
GAGE10_prt <- proteinToGenome(GAGE10_prt, edbX)


Running Snakemake in one single conda env

A good practice is to provide yml file along with your Snakefile. This file should contain all the enviroment definitions that are needed to run your pipeline. The minimum content of this file would be to install snakemake itself. Any other program used in your pipeline should be included as well. For reproducibility version numbers must be used.

So, if someone would like to use your pipeline, (s)he first creates a new enviroment based on your definition:

$ conda env create -f my_pipeline.yml

Afterwards activate this environment and run the pipeline:

$ conda activate my_pipeline
$ snakemake

The conda directive in your Snakefile is only neccessary if you need to use program which dependencies cannot be resolved within the main environment. If you start snakemake with the --use-conda parameter, snakemake will create a new environment based on the provided definition file and activate it for the rule.

  Especially since it is pulled down every time someone analyses data (sometime multiple times a day).

The environment creation is only done once. If you run the pipeline a second time snakemake can activate the environment diretcly without downloading and installing something again.

The tutorial I wrote some time ago on biostars about snakemake and conda might be useful for you.


Compare functional enrichment between genelists

To compare gene list based on gene ontologies you can use GOSemSim (Note, I'm a contributor to the package). You can select the semantic similarity to compare them, that is how the redundancies and hierarchical relations are counted. The main result you can obtain is similarity matrix for each comparison between gene lists.

However if you perform a GO term enrichment and then you compare the gene list of each GO term by looking at which GO do they belong it won't be much useful. I would suggest you to use BioCor (I'm the developer of it), which can use any pathway information (such as the ones of Reactome) to asses if they are related to the same pathways.


I was hoping to find some tool that already exists.

The short answer is, none exist yet. I see a number of hurdles to overcome first.1

For example, no one has yet characterized saturation for transcript discovery, i.e., the reads per cell needed to maximize transcript information recovery.  This is a typical stepping-stone as a sequencing approach matures, and until it's done, it makes it unlikely that most scRNA-seq datasets are going to be appropriate vis-à-vis statistical power, even if there was a tool.  That is, if people are only aiming for cell-type discrimination via gene-level markers, datasets will likely be underpowered to measure transcript level differences.

Another issue is the lack of any transcript-level scRNA-seq simulators. All the simulators out there that I know of produce gene-level matrices.  Having a transcript-level simulator is crucial for benchmarking and characterizing the performance of differential expression calling methodologies.  Without anything to test against, that leaves a huge blindspot in terms of estimating the uncertainty (false positive rates, sensitivity) that any tool might have.  However, since some scRNA-seq models (e.g., Velten, et al. 2015, or Lopez, et al., 2018) are generative, perhaps this is shouldn't be too difficult to adapt to generating FASTQs.

Simulators are fine and all, but ground truth is probably the key hurdle of the enterprise right now. While there are definitely groups using smFISH to distinguish isoforms quantitatively on a single-cell level, there doesn't appear to be enough smFISH data out there to measure a transcriptome-wide scRNA-seq method against. Hopefully this will change in the coming year or two, as results from more highly-multiplexed single-molecule technologies become available.

1 This is just my opinion formed as I try to keep up with the literature, so I could be missing things.  Comments/corrections certainly welcome!


How to visualize genome track of gene in specific cell-lines?

You can use screen.encodeproject.org
It collected a lot of data from multiple cell lines or tissues. And you can export tracks to UCSC browser or other browsers.


How can I match sequences with IDs using python dictionaries?

Biopython can index a fastafile and you can access the sequences by their ID in a dict-like manner.
All that is left for you to do is extract the IDs from your file. (more Biopython)


How to handle control samples in CLIP-Seq

I've not handled CLIP-seq before, but we have done iCLIP and in this we had controls, but no input. @DevonRyan is correct, these control samples are similar to IgG pulldowns in ChIP seq.

The way we handled this was to call our peaks on isolated samples. The original peak calling alog from the first CLIP seq paper doesn't require input. Also I believe that clipper doesn't need input (although it will take it). We then removed any peaks from the sample that overlapped a peak in the control.

But most of our analyses we didn't do on peaks - we did lots of things like looking at metagenes, or classes of genes, or particular types of exon etc, and when we did this, because we were averaging over many sites, we could generate the pattern for the control as well as for the input, and often normalise the one to the other.


Creating a new reference genome B by changing genome A with mismatches from BAM file without long reads

You should be able to use something like Pilon to convert from one genome to another, assuming the changes are all fairly small local changes (but that's almost never the case).

I find your resistance to long-read data interesting. It's important to consider the cost of bioinformatics as well. Sequencing a phage genome with long reads will take less than an hour, and probably produce a few reads that cover the entire genome. With existing short-read data, those full-length reads should be able to be fixed up within a day.

I spent about half a year trying to fix a eukaryotic genome assembly that was assembled only from short reads (using Pilon, and L-RNA-scaffolder), and I was only able to reduce the number of contigs down by about half (and even at that, I wasn't convinced that contigs were being properly joined). With long reads and another couple of weeks of computer work, the contigs dropped down to about a tenth of their original number.

If you can use a pipette and have about half an hour spare, you can do rapid long-read sequencing on the MinION.


How to dump genes from GenBank in GFF3 format?

Turns out you can just grab the GFF3 from the NCBI's FTP site!

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/003/195/GCA_000003195.3_Sorghum_bicolor_NCBIv3/

See:

https://www.ncbi.nlm.nih.gov/genome/doc/ftpfaq/#files

Or "Access the data" on the right here:

https://www.ncbi.nlm.nih.gov/assembly/GCF_000003195.3


Delete all 4 lines of a fastq read from a fastq file using read ID

Don't do it: Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ.  For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.

If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example FastqGeneralIterator and dumping the parsed result back to FASTQ in a 4-line-per-record form.


How can I use my Myheritage DNA results file for further analysis?

A csv file with RSID, chromosome, position and result is enough for what you want to do and these are the core columns of a VCF (which is just a TSV with some headers describing how it was made).
Given you have ~700k rows I suspect that yes there will be genetic health-related SNVs (single nucleotide variant). Disclaimer; these data could include info pertaining to your health and your family's health, I strongly advise that you speak to a genetic counsellor to understand these things.
Gnomad is a population database so generally speaking (plenty of exceptions of course) if there is a SNV in there more than a few times its probably not causing disease. https://www.ncbi.nlm.nih.gov/clinvar/ and https://omim.org/ are example of disease databases.
Here's a pandas example with clinvar (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/ you'll need to match your version).
import pandas as pd
df = pd.read_csv(
    'clinvar.vcf.gz',
    sep='\t',
    comment = '#',
    header=None,
    names = ['CHROM','POS','ID','REF','ALT','QUAL','FILTER','INFO']
)
df.head()

CHROM   POS ID  REF ALT QUAL    FILTER  INFO
0   1   930188  846933  G   A   .   .   ALLELEID=824438;CLNDISDB=MedGen:CN517202;CLNDN...
1   1   930203  972363  C   T   .   .   ALLELEID=959431;CLNDISDB=MedGen:CN517202;CLNDN...
2   1   930248  789256  G   A   .   .   AF_ESP=0.00347;AF_EXAC=0.00622;AF_TGP=0.00280;...
3   1   930275  969662  T   G   .   .   ALLELEID=959432;CLNDISDB=MedGen:CN517202;CLNDN...
4   1   930336  843786  G   A   .   .   ALLELEID=824439;CLNDISDB=MedGen:CN517202;CLNDN...

Then lets say one of your 700k SNVs is at chrom 1, position 930336 and you have an A
print(*df.query('(CHROM == 1) & (POS == 930336) & (ALT == "A")')['INFO'].str.split(';'))

Gives a list from which you could pull out CLNSIG which here is Uncertain_significance
['ALLELEID=824439', 'CLNDISDB=MedGen:CN517202', 'CLNDN=not_provided', 'CLNHGVS=NC_000001.11:g.930336G>A', 'CLNREVSTAT=criteria_provided,_single_submitter', 'CLNSIG=Uncertain_significance', 'CLNVC=single_nucleotide_variant', 'CLNVCSO=SO:0001483', 'GENEINFO=SAMD11:148398', 'MC=SO:0001583|missense_variant', 'ORIGIN=1']

So you could easily parse your 700k SNVs and get their respective clinical info this way. Please think about this carefully before you do it though, there are many implications including life insurance and mental health! I'm of the opinion that this sort of data is useful but MANY people disagree. Further, this advice is general in nature only, I accept no liabilty for what you or anyone else chooses to do with this publicly available data.


Is nanopolish worth it since faster polishing software is available?

Nanopolish is necessary if you want to get high-quality consensus. racon etc don't use signal data. They can't achieve high quality. pilon at times doesn't work well. Even if it worked perfectly, it wouldn't help repetitive regions inaccessible by short reads.


One option is to treat all proteins that are not explicitly known to form a biofilm as proteins that are incapable of doing so, but this will likely result in many false negatives.

You could take a look at one-class classification and positive-unlabeled (PU) learning. Both are techniques specifically designed for these kind of problems where only one of the two classes is known. (In this case: only biofilm forming proteins are known.)


How to filter blast results from blast run against virulence factor database

The answer is it depends. 40% pident with ~100% qcovs (assuming length of protein is not very small) tells you that the virulence factor (VF) and your prokka predicted protein are likely homologs, they belong to the same (or closely related) protein family. For some research studies, you might annotate your prokka proteins as 'putative' VFs.

On the other hand, in a clinical context, even 90% pident might not be enough evidence to say these predictions are real VFs. A single amino acid change alters the protein and therefore potentially the function. You could investigate bioinformtiaclly how the changes impact the protein structure and domains etc. Ideally validate it experimentally, e.g. for a new putative E. coli toxin, clone the protein and express it in a commensal E. coli strain and test the toxicity of the new strain.

In practise, think about what question you're asking and who you want to convince. Try numerous thresholds and see how robust your results are to threshold changes (4,000 VFs seems like a lot for a 'few' genomes). You could also try the VirulenceFinder database. When annotating antibiotic resistance genes, if you use 80% pident and qcovs thresholds with a minimum protein length then most researchers will accept your predicted annotations as 'conservative', maybe the VF community has similar arbitrary but pragmatic standards.


Rosetta ab initio run gives 0.00 RMSD for generated 2000 models? (run still in progress)

RMSD is a measure of similarity between two structures. The relevant RMSD in protein structure prediction is between a model and the native structure. If you don't have a native structure you cannot have an RMSD in this sense. Ignore the fact it is 0.0, there is possibly some error that means you are calculating the RMSD between a model and itself.

And that's the difficulty of protein structure prediction... you don't know whether you are right or not! The structures you have generated are your results, and you'll need to find ways to filter or use them that are not reliant on RMSD. You could use estimation of model accuracy methods or clustering, for example.


Alignment QC differences between HISAT2 and Qualimap

89038751 pairs times 2 is 178 million reads. That times the approximately 97% alignment rate is around 172 million reads, which is what qualimap is reporting.

A "non-unique" alignment is one that can map to multiple places (either equally well, or where all of the alignments are of similar enough quality). Such alignments have a single "primary" alignment and one or more "secondary" alignments, which is where the "secondary" metrics come from.

For the qualimap output, there are additional metrics describing how reads overlap genes (or not). Of note are the "ambiguous" alignments, which are those overlapping multiple genes. "No feature" metrics, then, refer to alignments not overlapping any gene (see the qualimap documentation for further details).

My guess is that you have such a high "intergenic" rate because you've either used the wrong GTF file or have told qualimap the wrong strand setting. For unusual organisms one expects many reads to arise from unannotated genes. This is not the case for common model organisms, such as mice and humans, where there are very few remaining commonly expressed unannotated genes.


Like gringer already commented, the cause of your error is nicely given by the java program. Your default max heap settings are probably too low for the program. If you have more memory on your machine you can give max heap a higher value. See this post and answers on Stack overflow.

You'll probably need to use the -Xmxn argument in java.


Your command is fine. It is working as expected. Did you intend to include -X 5 in your command? That restricts the number of spots to 5. If you drop it, you will download the entire set of 2563898 spots.


Is there a tool that can take a protein's amino acid sequence and would display it's locus on the genome?

I think the mistake here is thinking that you need to use the sequence in order to find the genome co-ordinates. If these proteins come from a well annotated speices, you can easily do this with BioMart.

Go the the ENSEMBL web site.
From the top menu, select BioMart
Under "-CHOSE DATABASE-" select Ensembl Genes, and under "-CHOOSE DATASET-" select your species (eg. Human Genes).
Select "Filters" from the lefthand menu and then expand the "GENE" section.
Select the "Input external references ID list [Max 500 advised]" and from the drop-down select "PDB IDs" or the correct "Uniprot KB" ID.
Paste your IDs into the box.
Select "Attributes" on the left-hand menu. Making sure you've the "Fetaures" radio box select, you can select Chromosome and either Gene Start/Gene End or Transcript start/Transcript end from the GENE section. These coordinates will also contain the untranslated regions unfortunately.  You may also want to play about with the various things under "PROTEIN DOMAINS AND FAMILIES", particurlarly under "Protein Features", although I'm not sure what they all mean.
When you are happy with the selection, select the big "Results" at the top to get the coords.

You can also get out Sequences for both the full transcript and the CDS and known germline and somatic variants that fall within the boundaries of the genes.


The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/

You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:

$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses chr# format and Ensembl uses #. For this, you can use the -A parameter from featureCounts. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:

## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv


ScaleData() from Seurat causes crash on RStudio Cloud

It seems RStudio server doesn't know it's limitations. Currently, Projects are limited to 1GB of RAM. Using the R package ulimit which mirrors the function of ulimit and allows implementation of memory limits in linux I was able to set my own limitations within the R session. After setting the limit to 1GB I see the ulimit limitation being tripped when invoking ScaleData(). Looks like a case of exceeding available system resources.


Analysis of differential transcript usage (DTU)

I will take the liberty of giving one possible answers to my own question – but I’m very interested in other answers.

One analysis type that such data enables is the analysis of transcript switches with predicted potential consequences.

I myself have recently developed such a tool called IsoformSwitchAnalyzeR. IsoformSwitchAnalyzeR enables statistical identification (via DRIMSeq) of isoform switches with predicted functional consequences. The consequences analyzed can be chosen from a long list which includes gain/loss of protein domains, signal peptides changes in NMD sensitivity etc. The R package also enables easy visualization of isoform switches along with their consequences and it directly supports the output of Cufflinks/Cuffdiff, RSEM, Salmon and Kallisto.

Apart from enabling identification of interesting examples switch identification also enables systematic analysis of what genes are affected. For inspiration I will recommend one of my own articles (describing results obtained with IsoformSwitchAnalyzeR) as well Hector et al’s recent paper (which is not using IsoformSwitchAnalyzeR). Of particular interest and finesse is Hector’s analysis of how isoform switches can disrupt protein-protein interactions.

Looking forward to hear more ideas.


Mutations in sequences relative to reference (histogram)

You can use AlignBuddy to process alignments in GenBank format, which allows you to keep the metadata. From the command line:

$: alignbuddy input.gb --generate_alignment mafft | \
   alignbuddy --extract_feature_sequences 'PROT_FEAT_1' | \
   alignbuddy --percent_id

If you want to do a large one-to-many comparison and loop through a bunch of features, I'd recommend importing AlignBuddy into a custom python script and calling the relevant functions directly. You can then collect the relevant data and build a histogram with matplotlib.


Disk space error while aligning reads using STAR

STAR versions 2.6.0b and 2.6.0a are unstable,  as written Issues performing variant calling with GATK you are using version 2.6.0b. You should switch to version 2.6.0c which is stable.  Also avoid to perform multiple alignments in the same directory.


How to merge back the processes by sample after splitting by chromosome in Nextflow?

Use the groupBy operator. If the content of channel xasxa is:
["/path/to/workdir/sample1_chr1.tsv", "/path/to/workdir/sample1_chr2.tsv", "/path/to/workdir/sample2_chr1.tsv", "/path/to/workdir/sample2_chr2.tsv
You could probably do:
xasxa.groupBy{ it -> it.fileName.toString().split("_")[0] }
The result would be a channel with:
[sample1:[/path/to/workdir/sample1_chr1.tsv, /path/to/workdir/sample1_chr2.tsv, sample2:[/path/to/workdir/sample2_chr1.tsv, /path/to/workdir/sample2_chr1.tsv]]
This breaks if you have an underscore in your samplename before _chr*. More robust would be to initialise all channels with a sample value, so that each channel is in the form of [samplename, chromosome]. Then all you need is a groupTuple on your channel.