Error processing run file with bcl2fastq - Unable to open '.../RunInfo.xml' file for reading I'm attempting to demux a nextseq run which I've successfully done several times in the past on other runs. I uploaded the run from my PC to the server (Ubuntu 18.04.4 LTS) without issue. Then I ran:
```bcl2fastq --runfolder-dir /data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0```

From which I got the error:

'/data/201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml' file for
reading```

However, the RunInfo.xml file is present and readable:
$ ls -l
total 96
drwxrwxrwx 10 user user 90112 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF
-rwxrwxr-x  1 user user 1929 Nov 11 07:37 SampleSheet.csv
$
$ ls -l 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
-rwxrwxr-x 1 user user 28510 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
$
$ head 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
<?xml version="1.0"?>
<RunInfo xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" Version="4">
  <Run Id="201109_NB551126_0128_AH5G5KBGXF" Number="128">
    <Flowcell>H5G5KBGXF</Flowcell>
    <Instrument>NB551126</Instrument>
    <Date>201109</Date>
    <Reads>
      <Read Number="1" NumCycles="75" IsIndexedRead="N" />
      <Read Number="2" NumCycles="6" IsIndexedRead="Y" />
    </Reads>

Has anyone had this issue or know of a fix? Thanks!
Update: contents of run folder
total 2452
-rwxrwxr-x 1 stewart stewart   977 Nov 11 07:59 '2020-11-09__14_05_58_94F12309 GetBestZ Camera 1 Focus Red top surface.snr.csv'
....
-rwxrwxr-x 1 stewart stewart   244 Nov 11 07:59 '2020-11-09__14_08_54_FE05AB9F SafeParabola called by CalculateSignalStd.snr.csv'
-rwxrwxr-x 1 stewart stewart   253 Nov 11 07:59 '2020-11-09__14_08_54_FE0620BC SafeParabola called by FindPeakUsingParabolaFit.snr.csv'
....
drwxrwxrwx 2 stewart stewart  4096 Nov 11 07:59  Config
drwxrwxrwx 3 stewart stewart  4096 Nov 11 07:59  Data
drwxrwxrwx 3 stewart stewart  4096 Nov 11 08:07  Images
drwxrwxrwx 2 stewart stewart 12288 Nov 11 08:07  InstrumentAnalyticsLogs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  InterOp
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Logs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Recipe
-rwxrwxr-x 1 stewart stewart    48 Nov 11 08:08  RTAComplete.txt
-rwxrwxr-x 1 stewart stewart  6200 Nov 11 08:08  RTAConfiguration.xml
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  RTALogs
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead1Complete.txt
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead2Complete.txt
-rwxrwxr-x 1 stewart stewart   925 Nov 11 08:08  RunCompletionStatus.xml
-rwxrwxr-x 1 stewart stewart 28510 Nov 11 08:08  RunInfo.xml
-rwxrwxr-x 1 stewart stewart 26802 Nov 11 08:08  RunParameters.xml
-rwxrwxr-x 1 stewart stewart  1929 Nov 11 08:08  SampleSheet.csv```

I repro the error.
ls /ml

ls: cannot access '/ml': No such file or directory

If run
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

I suspect it will work.
At present you are attempting to call a directory /data, which doesn't exit. If you type ls /etc that will work. The difference between ./data and /data in a file structure is alot, but on a command line it appears minimal.

Wait ... if that doesn't work .. because I suspect the calling directory is data then the following will,
bcl2fastq --runfolder-dir ./201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

OR
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./data/SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

It depends on the path where you are making the call. You could sort this out using $PWD or just echo $PWD then copy and paste the path name.


Disk space error while aligning reads using STAR Hi am trying to align RNA-seq reads using STAR through following commands,

STAR --runMode alignReads --runThreadN 26  --outSAMtype BAM  SortedByCoordinate Unsorted --sjdbOverhang 100  --genomeDir wheat --readFilesIn G3_cleaned_R1.fastq G3_cleaned_R2.fastq --sjdbGTFtagExonParentTranscript wheat.gff3
Jun 05 13:53:59 ..... started STAR run
Jun 05 13:53:59 ..... loading genome
Jun 05 14:07:12 ..... started mapping
Jun 05 18:42:47 ..... started sorting BAM

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin   does not agree with the actual size on disk: 509186583   0   45

EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 648593301   0   47

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 660412146   0   46
Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 470931698   0   44

Jun 05 18:42:47 ...... FATAL ERROR, exiting
EXITING because of FATAL ERROR: number of bytes expected from the BAM bin does not agree with the actual size on disk: 653155717   0   43

Jun 05 18:42:47 ...... FATAL ERROR, exiting
Segmentation fault (core dumped)

Although there is much free space on my system as shown below:

> df -Th
  Filesystem     Type      Size  Used Avail Use% Mounted on
  devtmpfs       devtmpfs   63G     0   63G   0% /dev
  tmpfs          tmpfs      63G  136K   63G   1% /dev/shm
  tmpfs          tmpfs      63G  1.9M   63G   1% /run
  tmpfs          tmpfs      63G     0   63G   0% /sys/fs/cgroup
  /dev/sda3      ext4      3.6T  1.4T  2.2T  39% /
  /dev/sda2      ext4      9.8G   57M  9.2G   1% /boot
  tmpfs          tmpfs      13G     0   13G   0% /run/user/480
  tmpfs          tmpfs      13G   16K   13G   1% /run/user/1000

How this problem can be solved?

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 change the UMAP use in the dimplot and feature plot I would like to know how to change the UMAP used in Dimplot and FeaturePlot from Seurat: how we can get the x-axis and the y-axis like UMAP-1 and UMAP-2 if I want to use UMAP-4 and UMAP-5. How to modify the code? It is easy to change the PC by using DimPlot(object = pbmc_small, dims = c(4, 5), reduction = "PCA") but if I changed to reduction = "UMAP", I got the error

"Error in Embeddings(object = object[[reduction]])[cells, dims] : subscript out of bounds Calls: DimPlot Execution halted".

Does anyone have any idea how to solve this problem?

Try this:

library(Seurat)
pbmc_small <- RunUMAP(pbmc_small, dims = 1:10, reduction = "pca", n.components = 5)

names(pbmc_small[["umap"]])
#> [1] "UMAP_1" "UMAP_2" "UMAP_3" "UMAP_4" "UMAP_5"

# plot
DimPlot(pbmc_small, reduction = "umap", dims = c(4, 5))

This should work with Seurat version 3.1.5

packageVersion("Seurat")
#> [1] '3.1.5'


Syntenic gene browser I found that there is a Syntenic gene browser at GEvo. Do anyone know where could I find a similar browser for a local installation?

The UCSC genome browser supports local mirrors. You could make a local mirror that has the "TransMap Alignments" cross-species gene track. From the track description page:

  These tracks contain cDNA and gene alignments produced by the TransMap cross-species alignment algorithm from other vertebrate species in the UCSC Genome Browser. For closer evolutionary distances, the alignments are created using syntenically filtered LASTZ or BLASTZ alignment chains, resulting in a prediction of the orthologous genes in human. For more distant organisms, reciprocal best alignments are used.

They provide a helpful tutorial on setting up a local installation here.


How can I find out the code DRIMSeq uses to make dmplotProportions objects? I'm using DRIMSeq to do some analysis but I'd like to make my own figures from the results ideally using what they've done as a base.

I'm trying to find how they make their ggplot objects in the function plotProportoins but when I use findMethods('plotProportoins') it ends up with a method dmPlotProportoins.

It says in the manual that it's a ggplot object but I can't seem to find how they make it. Any ideas?

I was curious how they have done it as well, so here expanded @Liopis comment :

wget http://bioconductor.org/packages/release/bioc/src/contrib/DRIMSeq_1.6.0.tar.gz
tar -zxf DRIMSeq_1.6.0.tar.gz
# if you need to find it, or use grep
# find DRIMSeq -name "dm_plotProportions.R"
less DRIMSeq/R/dm_plotProportions.R


Populations genetics and dynamics of bacteria on a Graph Disclaimer:
I'm a physicist and I'm fairly new to Bioinformatics therefore somethings below may not make sense.

My purpose:
I would like to simulate genetic evolution of bacterial populations, implementing interactions among different bacterial strains on a graph.
My starting point is to generate a graph with a certain topology.
Each vertex $V_i$ should contain a population of $N_i$ individuals for which I choose a certain allelic profile from an MLST database. Population numerosity and allelic profile is assigned by random choice.
$$\\\\$$
At starting time each vertex is populated with a "pure" strain of $N_i$ clones, identified with a randomly choosen allelic profile.
Dynamics assumptions
I will assume that:

Topology does not change with time

Vertex dynamic flows as a Moran process:
At time $t$ , one individual is chosen randomly to reproduce and one
individual is chosen to die. The same individual can be chosen to reproduce and
then die. Thus, an individual has either zero, one, or two descendants. Zero and
two with equal probability $p_0 = p_2 = (N_i -1 )/N_i^2$ , and one with probability $p_1 = 1 − 2 p_2$.

Migration from one Vertex to one another is possible for nearest neighborhood with some fixed probability depending on choosen topology. I will assume that when migration occurs $I$ new individuals enters into $V_i$ and $O$ individuals exits from $V_i$, with $I$ and $O$ depending on origin vertex numerosity. Migrants carry allelic frequencies according to population they came from. When migration occurs a new population is established: allelic profile frequencies changes according to that.

Migration and Moran process evolution does not happen simultaneously; This allows me to guarantee that each population performs at least one evolution.

Questions:
(1) Is there something fundamental that I'm missing?
(2) Moran model ensures that the population size remains constant (apart from migration effects obviously). Is this appropriate in order to simulate dynamic in stationary growing phase or should I introduce some demographic effects?
(3) Which could be a proper "timing" for migration and evolution steps?

I assume you are simulating a null distribution. Are you investigation recombination??
My main advice is to use population genetic terminology rather than geomometry to describe your simulation (e.g. vertex is a migration event between allopatric populations). Migration is investigated between populations investigated via $F_{ST}$. You appear to be assuming clonality, and population movement appear to follow a 'stepping stone model' (I might have misunderstand the question).
Mutation rate ... I'll have a guess at $10^5$ per site per year. Migration rate depends on the bacterial species being modeled, eg. MRSA spread in a very similiar fashion to that described in your model (along a moterway), in contrast soil bacteria don't adhere to allopatry at all (when I looked at them).
Fixing one parameter (population size) whilst investigation a second (migration) would be normal to assess the parameter space initially. Bacterial populations do undergo selective sweeps and using an appropriate growth model I would look at a standard population growth model (e.g. the one Bob May investigated in chaotic interactions within a differential equation). There are several population growth models and fixing migration whilst assessing the model would seem sensible. The parameter space in the combined model of migration, space and population size becomes complex.

Okay I get your model now. Its fine, the thing I'd be cautious about is the population density of the bacterial host because this would skew migration across the nodes of your graph. A quick literature review will demonstrate the importance,  double paeroto log normal distribution and its cousin (the name escapes me .... ) are sort of key words that might flag up. Bacteria are not my active thing so I'm rusty.


how to evaluate neural networks I have made some architectures of neural networks in order to classify proteins in three categories. I have calculated true positive, true negative, false negative, false positive, Mathews, sensitivity, specificity, accuracy for each category. My question is what is the right way to evaluate the architecture. Except for the c-index, could I calculate a cumulative sensitivity, accuracy etc?or should I evaluate separately for each category? Thank you!

Wow, I would just go for accuarcy as the singular most important measure, but be careful of sampling bias.

If the classification training had an unbiased sample* and you have an accuracy of >90%, in general ANN that wouldn't be great, but in bioinformatcs thats really good. If your training comprises sampling bias that gets tricky and you either need to accomodate that in your accuarcy ROC etc ... or else retrain to eliminate the sample bias.

You can account for all of them of course by wrapping everything into a single analysis.

*, by this I mean that each of the 3 classifications comprises 33% in each training set. If this is badly skewed, this requires consideration and gets complicated. If it is precisely 1/3 each then 'accurary' is the key stat.


How to install genometools with anaconda? I tied to install genometools with anaconda:

conda install -c bioconda genometools

Though, after running that successfully the command gt is not available in my conda environment.

  The program 'gt' is currently not installed. You can install it by
  typing: sudo apt install genometools

While conda tells me the program is installed. Are there multiple packages called genometools? What is the difference and can I install the package providing gt with anaconda at all? There is also a package called genometools-genometools, but it only has two downloads. So, I am not sure that is what I need.

genometools-genometools seem to be the current version. By installing the genometools package I received the same error like you. I guess it's because the binary isn't copy to bin during installation.

genometools-genometools seems to work fine. And it has 2k (2000) and not 2 :)


Removing repeated reads from nanopore 1D² reads This is summary of discussions and a question that was posted by Devon O'Rourke on Twitter

Following Albacore basecalls on a 1D² library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D² script. I'd like to use the reads generated from these scripts in a genome assembly.

After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D² directories - namely, those 1D reads that generated 1D² reads!

How do people approach the choice of selecting which directory to pull the unique read from if you're planning on assembling a genome? I'd think that the 1D² reads are of higher quality and be preferred, but perhaps there are 1D² reads which are lower in quality than one of the two 1D pairs?

Just looking for advice on how folks approach 1D² libraries, fastq outputs, and subsequent mapping. For what it's worth I'm hoping to use Minimap2 to proceed with alignments but am happy to hear others thoughts on alternatives.

Update

I tried the first solution proposed by @gringer, but my filtered_1d_reads.fq file contains same number of reads as raw_1d_reads.fq file.

I'm pretty sure the problem is the naming convention difference between the linear and paired reads. linear reads always follow a read_id of alphanumeric (8)-(4)-(4)-(4)-(12), while paired reads are (8)-(4)-(4)-(4)-(20)-(4)-(4)-(4)-(12). When I looked at it in more detail, the paired read "read_id" naming is just a concatenation of the two linear "read_id"'s that are partnered. So that (20) value is really the last (12) of the first linear read, and the first (8) of the second linear pair.

You can exclude the 1D² reads by generate a list of the ids from the 1D² reads, then filtering those reads out from the 1D fastq files. Here's one way to do that using some accessory scripts that I have written:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | awk '{print $2}' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

In order to deal with the concatenated read IDs, it's necessary to switch to a language like Perl that supports the use of backreferences for the splitting. The following perl command (replacing the first line in the above procedure) inserts a line break after the 36th character in the 2nd value, then prints that value:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt

So a full sequence of operations would look something like this:

$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq
$ cat 1d2_reads.fq filtered_1d_reads.fq > unique_reads.fq


How to use pyopenms.FeatureFindingMetabo() for peak picking in pyOpenMS? I found this documentation on pyOpenMS mass spec library. I want to try different peak picking algorithms. Unfortunately, the docs are not very detailed and don't show how to use the other peak picking algorithms and the codebase itself is also quite cryptic for me.
Any idea how to use, for example, the FeatureFindingMetabo algorithm insead of the centroided one?
from urllib.request import urlretrieve
# from urllib import urlretrieve  # use this code for Python 2.x
gh = "https://raw.githubusercontent.com/OpenMS/OpenMS/develop"
urlretrieve (gh +"/src/tests/topp/FeatureFinderCentroided_1_input.mzML", "feature_test.mzML")

from pyopenms import *

# Prepare data loading (save memory by only
# loading MS1 spectra into memory)
options = PeakFileOptions()
options.setMSLevels([1])
fh = MzMLFile()
fh.setOptions(options)

# Load data
input_map = MSExperiment()
fh.load("feature_test.mzML", input_map)
input_map.updateRanges()

ff = FeatureFinder()
ff.setLogType(LogType.CMD)

# Run the feature finder
name = "centroided"
features = FeatureMap()
seeds = FeatureMap()
params = FeatureFinder().getParameters(name)
ff.run(name, input_map, features, params, seeds)

features.setUniqueIds()
fh = FeatureXMLFile()
fh.store("output.featureXML", features)
print("Found", features.size(), "features")

This is a draft answer, if the OP still need needs I will extend it later.
Me too I found the documentation and examples of pyopenms very scarce. At the same time  OpenMS is huge, sometimes you find key information not in various other sources instead of the official documentation. This is especially true if it comes to metabolomics.
We built a Python OpenMS lipidomics LC MS/MS preprocessing pipeline, contained in this file:
https://github.com/saezlab/lipyd/blob/master/src/lipyd/msproc.py
The top class is the one called MSPreprocess, it implements the workflow, the other classes implement its steps.
Briefly, you start from profile mode mzML files, then

Peak picking with PeakPickerHiRes results centroided data
FeatureFindingMetabo assembles the traces in a metabolomics
compatible way
A MapAlignment algorithm alignes the traces across
scans
A FeatureGrouping algorithm alignes the features across multiple experiments

Our software mentioned above (lipyd) outputs data frames from the m/z, RT and intensity.
About OpenMS I can recommend to read the C++ code and docs. The Python classes and methods are automatically built minimal wrappers, everything works the same way as in C++.


Find rsIDs for GRCh37 SNPs and rsIDs for GRCh38 SNPs and compare for overlap I have a long list of variant IDs that I generated as a result of research done one GRCh37 genomes (e.g. 13_28025615_G_C_b37). I want to get their rsIDs and compare them to a newly-obtained list of GRCh38 rsIDs that have also been retrieved from similarly formatted variant IDs (e.g. chr22_50473581_G_A_b38). How would I go about doing this? I've tried installing biomaRt as a package but I can't quite make heads or tails of it and it's not clear to me that it would support both versions of genome builds.

it's not clear to me that it would support both versions of genome
  builds.

You can specify what version of the database you are querying.

I'd start with biomart on the ensembl website, see if you can get it to do what you want, then work on getting biomaRt to do what you want.


Is it possible to do alignment within Python? Check variants against reference? I'm currently adding a few SNPs randomly into a FASTA within python using BioPython. In the following example from BioPython, I add an SNV at location "5"

http://biopython.org/DIST/docs/api/Bio.Seq.MutableSeq-class.html

from Bio.Seq import MutableSeq
from Bio.Alphabet import generic_dna
my_seq = MutableSeq("ACTCGTCGTCG", generic_dna)
my_seq[5] = "A"
print(my_seq)
## outputs ACTCGACGTCG

In my case, I'm using the reference hg38 as the input FASTA to manipulate.

After randomly inserting these SNPs, I would like to find a region whereby there are no SNPs within (let's say) +/- 50Kb, or the next nearest thing.

Is it possible to align the FASTA with variants against the reference within python? Otherwise, is there a quick way to check the variant FASTA against the reference FASTA, finding a region with no variants? What would be an efficient way to execute this task?

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).


How can I run a command for multiple files? I have a bunch of .vcf files in a folder and I want to run command below on all of them but doing that one by one manually is really painful. I am seeking for a way in terminal to do that simultaneously. The output would be a .txt file.

grep -v "##" /file.vcf | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > /file.txt

Can you help please?

Assuming all your vcfs are in the same folder:

for file in *.vcf; do
    grep -v "##" ${file} | awk '{print $1"\t"$2"\t"$2"\t"$4"\t"$5"\t"$6"\t"$7"\t"$8"\t"$9"\t"$10"\t"$11}' > ${file}.txt
done


Any user friendly way to find rare mutations in whole genome raw? Is there any user friendly way to find rare mutations in the individual human whole genome sequencing raw data? (from Dante, 30x coverage).

To be more specific, I want to find mutations from this paper:

https://docs.google.com/document/d/1EkRMuD6J0-zyMY3MegKhz7Th0hzhzVXg6DWg2Hm6l2Y/edit?usp=sharing
(it's very short, less than one page).

But I'm confused. In their paper they just talk about genes (RAD21, B3GAT2, SMC3, SCN11A , SCN5A, SCN9A, SCN10A, SCN11A, TRPA1), but not mutations? Can we find that diseased-mutations which they talking about in genes, which their list? Or is there not enough data in the paper/study for this?

And if there is not enough data, which data I need to request from authors?

Or those genes have kind a "gold standard sequences" and if that sequence different from standard - there is "diseased" gene?
(I have only very basic genetics and bioinformatics knowledges)

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

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

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

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

People use tools like GATK to do this.

Hopefully that's enough information to get you started.


How to replace sequence identifiers in a fasta with OTU IDs from another file? I'm pretty new to Unix and bioinformatics and having a hard time accomplishing the following. I have one FASTA file with sequences and headers, and one OTU ID mapping file (.txt) with OTU IDs and their sequences. The sequence header in the .fasta file partially matches the sequence listed in the OTU ID file. What I want to do is replace the current headers in the FASTA file with just the OTU IDs from the .txt file. The files look like something like this:
OTU_ID.txt
denovo123 U1019.1010.2020_20384
denovo453 U1472.2992.2020_20443
denovo763 U1928.1898.2020_20343

.FASTA
>U1019.1010.2020_20384 0000-AHDJ-38347_223
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
>U1928.1898.2020_20343 0000-AHDJ-24247_242
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

and my desired output is that the .fasta file would look like this:
>denovo123
GAAATACCGAGAGAGAAATACCGAGAGAGAAAGAAATACCGAGAGAGAAATACCGAGAAATACCGAGAGAGAAATACCGAGAGA
>denovo763
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

Is it possible to do this in the command line?

You can do it with bash quiet easy. First create a backup of your fasta file with
cp otu.fasta otu_backup.fasta

Then you can run:
while IFS=" " read -r otu fastaId rest;
do
sed -E -i "s/>$fastaId .+/>$otu/" otu.fasta; done < ids.txt;
done < OTU_ID.txt

What does the code?
With IFS=" " you declare a whitespace as field seperator and by reading each line of your OTU_ID.txt, you have the two variables otu and fastaId which hold the otu name and the id of the Fasta Header.
Within the while loop, sed is looking each loop for the corresponding id within a RegEx and substitute it with the corresponding otu name. By using the -i option, you overwrite the existing file. This is why you need to create a backup.


What does "fetching by region is not available for SAM files" mean? I am used to gzip/biopython solutions when dealing with sequencing data, but now I wish to switch to more elegant pysam. So I looked at the manual, but ran into quite bizarre troubles with the first couple of lines using my bam file

import pysam
samfile = pysam.AlignmentFile("3_Tms_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tms_b3v08_scaf000159'):
    print(read)
samfile.close()

returns ValueError: fetching by region is not available for SAM files. Well, the file is bam. I tried to google the error but the only hits I found are the lines in the source code of pysam that check if the file is bam/cram or sam, so somehow pysam thinks that my bam is a sam. How can I convince it otherwise? I have also noticed that the manual is for python 2.7, that's maybe where the problem comes from...

That isn't actually a bam file as John Marshall figured out. I am keeping the rest of my answer since it could be useful to someone else, but the issue here was that you had a compressed (bgzipped) sam file and not an actual bam file and that's why you were getting that error. When I sorted your file in preparation for indexing it, I converted to a bam which is why the rest of my answer worked.

You don't have the index file for the bam file you're using. I used this script on the file you linked to (changing the name so that it corresponds to the right file and a contig in that file):

#!/usr/bin/env python3

import pysam
samfile = pysam.AlignmentFile("3_Tce_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
    print(read)
samfile.close()

The directory I ran it in had:

$ ls 3*
3_Tce_1_mapped.bam

And I got the error you described:

$ foo.py
Traceback (most recent call last):
  File "/home/terdon/scripts/foo.py", line 5, in <module>
    for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
  File "pysam/libcalignmentfile.pyx", line 1107, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetching by region is not available for SAM files

However, indexing the bam file fixed it:

$ samtools sort 3_Tce_1_mapped.bam > 3_Tce_1_mapped.sorted.bam
$ mv 3_Tce_1_mapped.sorted.bam 3_Tce_1_mapped.bam
$ samtools index 3_Tce_1_mapped.bam
$ ls 3*
3_Tce_1_mapped.bam  3_Tce_1_mapped.bam.bai

$ foo.py | wc
227    2724   16725

So just sort and index your files before attempting to seek in them. Which makes sense since the index's job is primarily to allow seeking.


Finding SNP/Indel Markers of Interspecific Variation This is a revised version of my original question, edited to include the portions already answered by @Maximilian Press.
Goal:
Find SNP and indel alleles which are fixed between but lack polymorphism within two closely related species (I'm looking for the BLUE region of the Venn diagram)
Desired Outputs:

All SNPs/Indels of Species A which are exclusive to and lack variation within Species A

All SNPs/Indels of Species B which are exclusive to and lack variation within Species B

Scratch Materials:

Species A reference genome

Species B reference genome

Illumina short reads for Species A

Illumina short reads for Species B

My current plan is to...

Map Species B to Species A:

minimap2 -ax asm5 referenceA.fasta referenceB.fasta > BtoA.sam

Find all SNP/Indel alleles of Species B relative to Species A (Represents the entire circle in the above Venn diagram):

bcftools mpileup -f referenceA.fasta -s BtoA.sam | bcftools call -o B_compared_to_A.vcf

Find all SNPs/Indel of Species B reads relative to Species B reference (Represents the red region of the Venn diagram):

bcftools mpileup -f referenceB.fasta -b bam_listB | bcftools call -o B_compared_to_B.vcf

Perform a liftover on "B_compared_to_B.vcf" to correct for coordinate differences between Species A and Species B references (I haven't completely figured out how to do this - I want to make sure I have the general pipeline nailed down before I get into the weeds)
Write a program which extracts all loci found in "B_compared_to_A.vcf" but NOT in "B_compared_to_B_liftedOver.vcf" (This should yield the BLUE region of the Venn Diagram
Repeat reciprocally to get alleles exclusive to Species A

Does the logic for this new pipeline check out?

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.


Single cell RNAseq cell cluster (true cluster or sub cluster) I am trying to run seurat on ~5000 sinngle cells. I am expecting minimum 15 cell types to be present in the data. I tried to runn it with multiple conditions;I can see there is 27 clusters; I believe many of them might be sub clusters. Is there a way to know which clusters are actually a sub-clusters.

As suggested by the members I have used clustree to check the cluster stability also changes the clustering parameter.

I wouldn't be surprised if the 15 cell types that you are expecting were mostly characterized via "protein-based" methods; antibody staining, flow cytometry using fluorophore coupled antibodies, ... The data that scRNA-seq is RNA based and should not be expected to provide the same information as protein based assays. Moreover, the scRNA-seq data, especially those from microfluidic systems like 10x, are sparse and might not be enough to resolve closely related cell types/states such as those of T cells.
On top of the aforementioned, evaluating clusters or clustering stability is a difficult task and in my experience, quite some cluster evaluation metrics cannot be easily applied to single cell data simply because of the shear size of the data. Having asked your question and looked for answers (and I am still doing that), I am doing the following:
i) Try to use transcriptomic markers as much as possible. Even then a lack of a marker would not mean much in terms of scRNA-seq, it might very well not be detected simply due to sensitivity
ii) Use clustree do select the most plausible number of clusters (resolution parameter in Seurat). The package is compatible with Seurat and some other scRNA-seq packages.
iii) Use silhouette width as a metric for clustering. This is computationally expensive, however, should be fine for your 5000 cells. For larger cell numbers, I use another package that approximates this (will add a link but first I have go through my scripts for the name).
iv) Check if batch effects are accounting for what you call "sub-clusters". For example clusters 0, 14, 24 and 25 above might correspond to different samples processed in different days.


Interpreting this plot from GSEA I have RNA-seq for two groups of patients: Responders to chemotherapy (n=9) versus non-responders to chemotherapy (n=24)
I have done DESeq2 and obtain fold change of responders versus non-responders in which positive fold change means down-regulation in responders (if I am not wrong)
> head(df)
    gene       FC
1 286499 7.116071
2 260436 6.691189
3   2267 5.751322
4   5968 5.522079
5   4103 5.328008
6 442709 5.036087
>

I have done GSEA in pre-ranked way
got this plot for KRAS pathway

I interpret this as; KRAS pathway in being down regulated in responders
Am I right?
Please correct me
When I separated genes for up and down regulated
For up regulated genes I did not get any significant GSEA but for down regulated genes in responders I got this

SO does this mean that KRAS pathway is being up regulated in responders because I have used down regulated genes in responders?

In Deseq2, it depends on contrast eg. if your contrast <- c("Condition", firstC, SecondC) is then -ve is downregulated in FirstC and same for +ve. It is always better to take look at count data for the confirmation (responders versus non-responders samples).
[for more detail you can refer here][1]
[1]https://github.com/amarinderthind/RNA-seq-tutorial-for-gene-differential-expression-analysis/blob/master/RNA-Seq-DGE.R


How to find Genotype, from list of mutations:lineage list? I'm very new to python and I'm getting along nicely (I'd like to believe); however, I must be missing something here. I'm looking to read each and every file and compare each line against one another.
The aim of this is to see if any mutations are shared in any of the strains - to determine the genotype.
The inner-workings of the files are single spaced headers, as those you would find on GISAID.
os.chdir("C:)
with open('A D614G.txt', 'r') as file2:
with open('A E780Q.txt', 'r') as file3:
with open('A G476S.txt', 'r') as file4:
with open('A L18F.txt', 'r') as file5:
with open('A N439K.txt', 'r') as file6:
with open('A S477N.txt', 'r') as file7:
with open('A T478I.txt', 'r') as file8:
with open('A V483A.txt', 'r') as file9:

files = [file1, file2, file3,file4, file5, file6, file7, file8, file9]

it = itertools.permutations((files), len(files))

for x in it:
    print(x)

os.chdir("C:)
with open('tester.txt', 'wt') as file_out:
    for x in it:
        file_out.write()```

1 I recieve the following error on a very long repeated stretch of text

```<_io.TextIOWrapper name='A G476S.txt' mode='r' encoding='cp1252'>```

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)


How to get the uniparc ids for all members of uniref cluster in batch request? I currently use the 'UniProt website REST API' (example) to collect the UniParc IDs of all the members of the given cluster. Fetching each entry individually is very slow though, so I wanted to ask if there is a way that I can do this request in batches.
What I have already tried:
There are pages mentioning a batch retrieval of entries. As I understood it, all the examples only adress the ID mapping service. Here I have tried adjusting the given example to the following
import urllib.parse
import urllib.request

url = 'https://www.uniprot.org/uploadlists/'

params = {
'from': 'NF100',
'to': 'UPARC',
'format': 'tab',
'query': 'UniRef100_Q8RKD8'
}

data = urllib.parse.urlencode(params)
data = data.encode('utf-8')
req = urllib.request.Request(url, data)
with urllib.request.urlopen(req) as f:
   response = f.read()
print(response.decode('utf-8'))

which resulted in the output
From    To
UniRef100_Q8RKD8        null

I also checked the interface of this service on the website and there I can't seem to map NF100 to UPARC, so I'm assuming that is why I'm not getting a result via the API.
I also know that the UniRef Id contains either a UniParc ID or the UniProtKB accession, so I could use the UniParc ID directly or fetch it for the UniProtKB via the proteins api with the uniparc endpoint, but from what I can see I only get the UniParc ID/entry of one of the members.

I found a way to batch the mapping from uniref cluster ids to the ids of all members. The url
https://www.uniprot.org/uniref/?query=UniRef100_A0A003&format=tab&columns=id,members
can be used to get the link from one uniref Cluster to it's member ids. Multiple uniref cluster can be fetched at once by connecting the ids with a +OR+ like so https://www.uniprot.org/uniref/?query=UniRef100_P99999+OR+UniRef100_P31946&format=tab&columns=id,members
Because the uniref cluster ids are not stable I ended up using the following syntax, where instead of querying the UniRef cluster by it's id, I query by one of it's members, namely the one that was used to create the current uniref cluster id: https://www.uniprot.org/uniref/?query=*&format=tab&limit=10&columns=id,members&fil=identity:1.0+AND+member:(P99999+OR+P31946)
In the notes section of the 'Retrieve / ID mapping' documentation of the uniprot they talk about fetching chunks up to 20.000 identifiers at a time. I am not sure if I use the service they are talking about though, so I'm keeping my requests at 100 identifiers per request and that seems to work.
One thing I noticed was that fetching results via the uniref endpoint is faster for small sets of identifiers (<20) since the fetch time of the batch request is similar for smaller and larger (<100) sets. So batching only makes sense for more than 20 identifiers at a time. (Disclaimer: I have only done very basic tests on this)
Once I have all the members of a cluster I can use the proteins/uniparc endpoint of the proteins api to fetch the uniparc data. For uniparc accessions (starting with the string 'UPI') I use the GET /uniparc/upi/{upi} endpoint. For the other accession, which should all be Uniprot accessions, I use the GET /uniparc/accession/{accession} endpoint.
This doesn't answer the question of a direct mapping between the uniref ids and uniparc member ids, but it is a workable solution.


Sleuth: transcripts with beta close to 0 are considered differentially expressed in a likelihood-ratio test I'm comparing the results that I obtain when doing a DE analysis with the Wald test and the likelihood-ratio test. One the thing that I've noticed is that there are many genes with 'beta' close to zero that are considered differentially expressed between the conditions.

I know that the likelihood-ratio test does't use the beta values to calculate the p-values, but I find it strange that transcripts with similar expression values between the conditions are being considered differentially expressed.

Volcano plot (Wald test q-values):

Volcano plot (likelihood-ratio test q-values):


how could i find uniprot ac in blast? I have made a search in p-blast.
I would like to search for the protein results in uniprot. how can I do it directly?
I can go to GENBANK but i would like to go to UNIPROT.
Thank you!!!

UniProt also provides a BLAST service (https://www.uniprot.org/blast), or you could use the Expasy BLAST too, at https://web.expasy.org/blast.
If you have performed your similarity search at NCBI or elsewhere, you can map the identifiers from that result to UniProt by using the UniProt IDmapping service at https://www.uniprot.org/uploadlists. This can also be used programmatically, cf https://www.uniprot.org/help/api_idmapping


Customizing colors in geom_arc_bar (ggforce) I am trying to plot scatter plot of pie charts using ggplot2 and ggforce.
Following is my code

library(data.table)
library(ggplot2)
library(ggforce)
library(dplyr)

df <- fread("data_frame_pie_plot.csv", header = TRUE,
        data.table = FALSE, stringsAsFactors = FALSE)

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

plot_obtained
Link_to_data Data_frame.
Can I fill the colours manually (discrete(based on discrete_fill column)
and continous (based on conti_fill column))?
I have tried the following codes.

##================================================================================##
# Code 1
my_colors <- c("red","white","blue")
ggplot(df, fill = sample_mut_colour) + req_data +
  theme_no_axes() + scale_color_manual(values=my_colors)

##================================================================================##
# Code 2
my_colors_1 <- c("red", "white")
my_colors_2 <- c("blue", "white")

req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill,
               color = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 3
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = case_when(discrete_fill == 1 ~ my_colors_1,
                                 discrete_fill == 2 ~ my_colors_2)), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##================================================================================##
# Code 4
req_data <- apply(df, 1, function(x){
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount,
               fill = discrete_fill),
           palette = case_when(discrete_fill == 1 ~ "red",
                               discrete_fill == 2 ~ "black",
                               TRUE ~ "white"), stat = "pie")})

ggplot(df) + req_data + theme_no_axes()

##=================================================================================##

But they do not seem to work. Any idea what is wrong in the code?

For a discrete scale your Code 1 example was nearly there.  Just use scale_fill_manual instead of scale_color_manual.  For more control over how the discrete values map to color values you can name the vector, too.

my_colors <- c("red","white","blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)
# or, if you want consistent output independent of the order of
# df$discrete_fill:
my_colors <- c("1" = "red", "2" = "white", "dummy" = "blue")
ggplot(df) + req_data +
    theme_no_axes() + scale_fill_manual(values=my_colors)

You can also just use a single call to geom_arc_bar for the whole data frame (rather than a list of call results from each row) and you'll get the same effect.

ggplot(df) +
  geom_arc_bar(aes(x0 = x0, y0 = y0, r0 = 0, r = 0.45,
               amount = amount, fill = discrete_fill), stat = "pie") +
    theme_no_axes() + scale_fill_manual(values=my_colors)


Meaning of BWA-MEM MAPQ scores Does anyone know what the MAPQ values produced by BWA-MEM mean?

I'm looking for something similar to what Keith Bradnam discovered for Tophat v 1.4.1, where he realized that:

  0 = maps to 5 or more locations

  1 = maps to 3-4 locations

  3 = maps to
  2 locations

  255 = unique mapping

I'm familiar with the notion that MAPQ should be theoretically be related to the probability of an "incorrect" alignment (10^(-MAPQ/10)), although this is vaguely-specified enough that aligners in practice tend to actually just use something like the above.

Note that this is different than the BWA MAPQ scoring interpretation, because BWA-MEM gives MAPQ scores in the range $[0,60]$, rather than $[0,37]$ as has been established for BWA.

First of all, if you want to understand mapping quality (mapQ), ignore RNA-seq mappers. They often produce misleading mapQ because mapQ is not important to RNA-seq anyway.

Strictly speaking, you have two questions, one in the title: the meaning of mapQ; and the other in a comment: how mapQ is computed. On the meaning, mapQ is nearly the same as baseQ – the phred scaled probability of the alignment/base being wrong. It often amuses me that we question mapQ but take baseQ for granted. BaseQ is also scaled and discretized differently; even fewer people know how Illumina/pacbio/nanopore/historical sequencers estimates baseQ.

On the second question, Section 2 of the MAQ supplementary explains the theoretical aspects of mapQ, which is still correct today. Briefly, mapping quality consists of three components: 1) the probability of contamination; 2) the effect of mapping heuristics and 3) the error due to the repetitiveness of the reference. Only 3) can be modeled theoretically.

In case of bwa-mem, if we assume the matching score is 1, type-3 error is estimated with:

$$10/\log10\cdot[\log4\cdot(S_1-S_2)-\log n_{\rm sub}]$$

where $S_1$ is the best alignment score, $S_2$ is the second best and $n_{\rm sub}$ is the number of suboptimal alignments. Factor $\log 4$ comes from the scoring matrix. Factor $10/\log10$ is the Phred scale. This equation assumes gap-free alignment and is very close to Section 2.5.2 in the MAQ supplementary. It is ok-ish for short reads, but often overestimates for long reads. I am not aware of a practical approach in general cases. In addition to this method, you can estimate mapQ by read simulation: just try to find a function such that it fits empirical mapQ. Some have tried machine learning, too.


How are different forcefields developed/evaluated? The Rosetta software can use a variety of different energy functions:

Ref2015
GenPot
Beta-Nov16

How are these energy functions developed and validated? Is it common to do a single protein design with several of these energy functions?

Some force-fields rely solely on theoretical physical principles, while others are tweaked specifically for the goal of accurately depicting certain structures like sheets and helices or movement/diffusion aspects at different theoretical temperatures of a set of well-studied proteins. (Top-Down vs. Bottom-Up)
I haven't done MD myself but I am in an MD group and people seem to normally stick to a single force-field for a given MD project but it happens that they do very short simulation runs with force-fields just as a test.
They definitely don't plan on doing the full simulation on different force-fields unless the comparison of those is their main goal.


Adding new columns to this data frame I have a dataframe like this

> dim(anno_maf)
[1] 32786   129
> head(anno_maf[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Start_Position
1:                    a       chr1       27105930
2:                    a       chr1        3342771
   End_Position
1:     27105930
2:      3342771
>

I have a list of genes like

gene name
PTEN
PTEN
MAP3K1
MALAT1
MED12
RUNX1
PIK3CA
FBXW7
ARID1A
MED12
TAB1

I want to add a column in my dataframe named COSMIC, if a gene in gene list exists in my dataframe Hugo_Symbol column, put a yes for this gene if not put a no; Finally my data frame has an extra column as COSMIC in which I have yes or no in Hugo_Symbol column based on the existence of genes in my gene list

Can you help me please?

Thank you

You can use ifelse():

anno_maf$COSMIC <- ifelse(anno_maf$Hugo_Symbol %in% your_gene_list,
                          yes = "yes",
                          no = "no")


Is there a point in recalibration of scores for variant calling? The most variant calling pipeline GATK include a Base Quality Score Recalibration (BQSR) which requires a list of known variants. Recently, some work has been done for reference-free recalibration of scores as well: Lacer and atlas, which is motivated by making the most for aDNA and low coverage datasets.

The importance for aDNA is explained in this lecture, but it is not clear to me if / how is important BQSR is for fresh DNA samples with decent (>15x) coverage. Especially when I work with non-model organisms and I can not simply use the standard tools.

How big an impact does recalibration of scores have on variant calling? Is there a rule of thumb for which it is / it is not worth the effort?

I personally don't think BQSR has a huge impact on variant calling, but you don't really need to guess. If you run GATK BQSR, it outputs a table and charts of exactly how much quality scores are adjusted. The adjustment will vary depending on the position in the read and genomic context (previous and following base). In my experience, the difference is a few points at most, but it's certainly noticeable.

GATK recommends BQSR for both genome and exome data, which is normally much higher than 15x.


Median string problem & multiple sequence alignment I read about the median string problem as an introduction to the multiple sequence alignment, however none of the MSA algorithms used seems to be using the idea of finding the median string.
To my understanding from the consensus (median) string, it would be possible to retrieve the actual alignment by simply doing the pairwise alignment of the consensus string with each of the strings.
Why is this approach not used? Are there any algorithms that use it?
And related to this, what kind of tasks in bioinformatics/biology are median/consensus string algorithms used for?

Why is this approach not used? Are there any algorithms that use it?

As shown in this paper, the MEDIAN STRING problem is NP-Complete on alphabets of carnality greater than 3. While the problem isn't used directly in approaches due to the complexity, the idea of trying to create a consensus sequence between pairs of input and using that as building blocks is common. These are considered "progressive alignment methods".


Convert PAF format to SAM/BAM format I have a bunch of PAF files resulting from the alignments of fastq files on a reference genome with minimap2. I would like to convert them into SAM/BAM format so I can use samtools mpileup to recover the read coverage at filtered positions

So far I've looked at paftools.js and found sam2paf but no paf2sam

I am open to suggestions but I'd like to avoid mapping the reads again.

Edit

I have opened an issue on minimap2's github

Answer from Heng Li :

  You need to rerun alignment. No tools exist to convert paf to sam.

why there is no paf2sam in paftools.js

  It is technically possible, but someone needs to write the code.

The paf format does not have the CIGAR string and base quality information that are typically present in a BAM/SAM file. While it is in principle possible to write a paf2sam converter, it's going to be very difficult to regenerate the CIGAR string with all the indel information (and maybe it's probably impossible to do so) that you'll need in a BAM/SAM file. If your motivation is to get coverage information, it might make more sense to convert the paf file into a bed file, before using a tool like bedtools to convert it into a coverage track.


How to indicate the END of a haplotype block in VCF? In VCF I know how to indicate that two genotypes are in the same phase by using consecutive "0|1" and "1|0" genotype fields, for example. However, how do I deal with the case that the first two genotypes are in phase and the second 2 are in phase, but there is no phase between the 2nd and 3rd?

According to the format I think I need to write:

chr1 100 . A T . . . GT 1|0
chr1 200 . A T . . . GT 0|1
chr1 300 . A T . . . GT 1|0
chr1 400 . A T . . . GT 0|1

but that implies that the SNP at 200 and the SNP at 300 are in phase, when in fact they are not.

How to indicate the END of a haplotype block in VCF?

According to the VCF specification, haplotype blocks are defined by the "Phase Set" [PS] tag:

  PS : phase set. A phase set is defined as a set of phased genotypes to which this genotype belongs. Phased genotypes for an individual that are on the same chromosome and have the same PS value are in the same phased set. A phase set specifies multi-marker haplotypes for the phased genotypes in the set. All phased genotypes that do not contain a PS subfield are assumed to belong to the same phased set. If the genotype in the GT field is unphased, the corresponding PS field is ignored. The recommended convention is to use the position of the first variant in the set as the PS identifier (although this is not required). (Non-negative 32-bit Integer)

A good example of this can be found on a 10X Genomics support page about phasing:

CHROM   Pos REF ALT GT  PS
chr1    1000    A   C   0|1 1000
chr1    1010    T   G   1|0 1000
chr1    2000    C   T   0|1 2000
chr1    2005    T   G   0/1 2000
chr1    2008    G   C   0|1 2000

  In this example we have two phase blocks, denoted by PS=1000 and PS=2000. PS=1000 spans position 1000-1010, and PS=2000 spans position 2000-2008. In PS=1000, haplotype 1 contains the REF A allele at position 1000, and the ALT G allele as position 1010, while haplotype 2 contains the ALT C allele at position 1000 and the REF T allele at position 1010.

Note that this phase set identifier needs to be a single numeric value, which means that sub-blocks and overlapping blocks cannot be defined by this field.


How to get rows with similar values in two different columns using command line? I have following example data from a large file.
cont100 1128 1125
cont1005 3642 3642
cont1006 103 19
cont1037 3146 3146
cont104 895 890
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577
cont1078 43 42
cont1081 727 2

I need to get only the rows with similar values in column 2 and column 3. The output should be like this.
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Also I need to get the inverse of that. That means I need to get the rows with values that are not similar in the columns. The output should be like this.
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

I think it is easy to do with awk, e.g. see here.
For your example, with your data in file:
%  awk '$2 != $3 {print $0}' file
# outputs
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

%  awk '$2 == $3 {print $0}' file
# outputs:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Note that where you say "similar" I am interpreting that to mean "exactly the same", because that is the case in your example.


Error while performing "lift-over" of genomic coordinates between two species using Ensembl Compara API I'm trying to convert coordinates between the reference mouse genome and CAST/EiJ mice strain. Since Ensebml already has alignment available in the browser, I decided to use Ensembl perl API to lift over the alignment for the regions I'm interested in. The problem is that it doesn't seem to be working. I'm getting following error message:

-------------------- EXCEPTION --------------------
MSG: The database do not contain any CACTUS_HAL data for Bio::EnsEMBL::Compara::GenomeDB=HASH(0x557db0d661e8)!
STACK toplevel coord.pl:42
Date (localtime)    = Mon Oct 15 07:42:17 2018
Ensembl API version = 94
---------------------------------------------------

which I assume means that I'm either searching for the wrong alignment record or the alignment itself is not available. Following is the perl script which I modified from the Ensembl Compara API tutorial:

use strict;
use Bio::EnsEMBL::Registry;
use Bio::EnsEMBL::Utils::Exception qw(throw);
use Bio::AlignIO;
use Getopt::Long;

my $species = "murinae";
my $targetSpecies = "mus_musculus_casteij";
my $alignment_type = "CACTUS_HAL";
my $output_file = undef;
my $coord_system = "chromosome";
my $seq_region = "14";
my $seq_region_start = 5000000;
my $seq_region_end = 5010000;
my $output_format = "clustalw";
my $genome_dbs;

use Bio::EnsEMBL::Registry;
Bio::EnsEMBL::Registry->load_registry_from_db(
    -host => 'ensembldb.ensembl.org',
    -user => 'anonymous',
    -port => 5306);

my $genome_db_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomeDB');

throw("Cannot connect to Compara") if (!$genome_db_adaptor);

my $this_genome_db = $genome_db_adaptor->fetch_by_name_assembly($species);
my $target_genome_db = $genome_db_adaptor->fetch_by_name_assembly($targetSpecies);
push(@$genome_dbs, $this_genome_db);
push(@$genome_dbs, $target_genome_db);

my $method_link_species_set_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'MethodLinkSpeciesSet');

my $method_link_species_set =
$method_link_species_set_adaptor->fetch_by_method_link_type_GenomeDBs(
      $alignment_type,
  $genome_dbs);

throw("The database do not contain any $alignment_type data for $target_genome_db!")
    if (!$method_link_species_set);

my $slice_adaptor = Bio::EnsEMBL::Registry->get_adaptor($species, 'core', 'Slice');

throw("Registry configuration file has no data for connecting to <$species>")
if (!$slice_adaptor);

my $query_slice = $slice_adaptor->fetch_by_region('toplevel', $seq_region, $seq_region_start, $seq_region_end);

throw("No Slice can be created with coordinates $seq_region:$seq_region_start-".
    "$seq_region_end") if (!$query_slice);

# Fetching all the GenomicAlignBlock corresponding to this Slice:
my $genomic_align_block_adaptor = Bio::EnsEMBL::Registry->get_adaptor(
    'Multi', 'compara', 'GenomicAlignBlock');

my $genomic_align_blocks =
$genomic_align_block_adaptor->fetch_all_by_MethodLinkSpeciesSet_Slice(
      $method_link_species_set,
  $query_slice);

my $all_aligns;

# Get a Bio::SimpleAlign object from every GenomicAlignBlock
foreach my $this_genomic_align_block (@$genomic_align_blocks) {
    my $simple_align = $this_genomic_align_block->get_SimpleAlign;
    push(@$all_aligns, $simple_align);
}

# print all the genomic alignments using a Bio::AlignIO object
my $alignIO = Bio::AlignIO->newFh(
-interleaved => 0,
-fh => \*STDOUT,
-format => $output_format,
    -idlength => 10
);

foreach my $this_align (@$all_aligns) {
    print $alignIO $this_align;
}

exit;

Any pointers on what I'm doing wrong ?

While I wouldn't consider this as the best answer or solution - I decided to use Ensembl REST api to extract the coordinates of alignment blocks.

The quick fix - in python looks like this:

#This is the main function which uses Ensembl REST API to access the Enesembl Compara database
#And retrive the alignment blocks which correspond to the MM10 standard.
#The function accepts chromosome, start and array and returns the alignment blocks
#reported by the API. Since the blocks can be on different chromosomes/and/or/contigs, it
#returns the whole output in form of dataframe - containing chromosome, start and stop of the
#alignment blocks corresponding to input and all the other alignment blocks corresponding the
#16 mice species, if available.
#!!!!!!!----------Caution start----------!!!!!!!
#The API is a bit buggy, so it is advised not to request for specific species.
#!!!!!!!-----------Caution end-----------!!!!!!!
def liftover(chrom, start, end):
    mmStart, mmEnd, mmchr = [], [], []
    caStart, caEnd, cachr = [], [], []
    #ryStart, ryEnd, rychr = [], [], []
    server = "https://rest.ensembl.org"
    ext = "/alignment/region/mus_musculus/" + str(chrom) + ":" + str(start) + "-" + str(end) + ":1?species_set_group=murinae;method=CACTUS_HAL;content-type=application/json"
    headers={ "Content-Type" : "application/json"}
    r = requests.get(server+ext, headers={ "Content-Type" : "application/json"})
    if(r.ok):
        json_string = json.dumps(r.json())
        jData = json.loads(r.content)
        results = r.json()
        print("The response contains {0} properties".format(len(jData)))
        print("\n")
    else:
        r.raise_for_status()

    for element in jData:
        element.pop('tree', None)

    for element in jData:
        for key, value in element.items():
            for z in value:
                if (z["species"]=="mus_musculus"):
                    mmStart.append(z["start"])
                    mmEnd.append(z["end"])
                    mmchr.append(z["seq_region"])
                if (z["species"]=="mus_musculus_casteij"):
                    caStart.append(z["start"])
                    caEnd.append(z["end"])
                    cachr.append(z["seq_region"])
    if (len(mmStart) != len(caStart)):
        print("Length mismatch between species blocks - proceed with caution!")
    return [min(mmStart), max(mmEnd), min(caStart), max(caEnd), len(mmStart)]

Again, it depends on how alignment blocks are reported - and requires further digging to resolve which mm10 alignment block resembles CAST alignment block.


Can a customized GRCh38 .gtf file be used with any of the GRCh38 released patches? I got a customized GRCh38.79 .gtf file (modified to have no MT genes) and I need to create a reference genome out of it (for 10xGenomics CellRanger pipeline). I suspect that the .79 part is the Ensembl number, which according this ensembl archive list is paired with the GRCh38.p2 patch.

Should I use this patch's fasta file, or it would be fine to use any of the GRCh38 patches?

You can see the exact assembly by looking in the README file in Ensembl's genome download page. ftp://ftp.ensembl.org/pub/release-89/fasta/homo_sapiens/dna/README

As you can see the current assembly is GCA_000001405.25. Ensembl 79 uses version GCA_000001405.17. You can expect some differences between assembly versions.

If you take a look at https://www.ncbi.nlm.nih.gov/assembly?term=GRCh38&cmd=DetailsSearch

You can see that 20,466,394 bp have been added since 17 -> 25. This is only about 0.6% change so I would guess there is not that much difference between these versions. However, if I were you I would use the corresponding version of Ensembl. Even better would be using the latest version of Ensembl, as not only the genomes primary sequence but also the positions of genes can change between versions (https://doi.org/10.1093/bib/bbw017).


Error in PSI-BLAST when running psipred locally I am very new to bioinformatics and Linux and am trying to run psipred locally on Ubuntu 18.04 but cannot get PSI-BLAST to work. I have downloaded BLAST legacy version 2.2.26 and unpacked it, uniref90 database and psipred version 4.01

I installed psipred according to instructions:

cd psipred/src
make
make install

And this seems to work fine

I used the following code to create a database, which proceeded without errors:

makeblastdb -in uniref90.fasta -dbtype prot -out uniref90

In runpsipred I changed the directories:

# The name of the BLAST data bank
set dbname = /home/hilda/uniref90

# Where the NCBI programs have been installed
set ncbidir = /home/hilda/blast/bin

# Where the PSIPRED V4 programs have been installed
set execdir = ./bin

# Where the PSIPRED V4 data files have been installed
set datadir = ./data

I then try to run psipred:

tsch ./runpsipred example/example.fasta

And get the error

Running PSI-BLAST with sequence example/example.fasta ...
FATAL: Error whilst running blastpgp - script terminated!

I am not sure how to solve this as I do not understand what is going wrong. As I said I am new to Linux. I use the command tsch to run psipred (otherwise I get the error "command not found"), but not to install the database, could this be a problem? Or have I missed something in the setup process? Running psipred does create a copy of the example.fasta in the psipred folder named psitmp4286007f010.fasta, but that is all that seems to happen.

I would greatly appreciate any suggestions on how to make it work!

Edit:
I tried to run with BLAST+ (which seemed to work without tsch) and get the following error:

BLAST Database error: No alias or index file found for protein database [/home/hilda/uniref90] in search path [/home/hilda/psipred::]

I have also tried setting dbname to home/hilda/uniref90/uniref90 and got the same results. there is a uniref90.pal file in the uniref90 directory, which I understand is the alias file? Other forum answers seem to use pfilt to index (?) the database, but there is no pfilt in my psipred directory. Any ideas how to fix the problem?

This error could be due to the fact that you are using legacy blast, have you tried using BLAST+ instead?
If you wish to use PSIPRED with BLAST+, then use the runpsipredplus script
here, rather than the normal runpsipred script.


Why do I obtain very strange results with DESeq2? I am using DESeq2 to perform a differential expression analysis, but I obtained very strange results for some genes. For some genes, I have very high log2FoldChange with very low p-value as displayed in the following table.

When I look the normalized expression values for these genes.  I can see that these genes are expressed in only one single sample (condition1= samples from 1 to 7, condition2 = samples 8 to 13).

Could you please tell me how to avoid DESeq2 to detect these genes as differentially expressed ?

I am using to following commands:

dds                      <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable, directory = directory, design= ~ condition)
dds                      <- DESeq(dds)
res                      <- results(dds, contrast=c("condition",condition1,condition2))

I'm guessing here, but I suspect this is at least in part due to the fact that all your samples in condition 1 are not expressed at all.

Formally you have an infinitely large fold-change between condition1 and condition2. Most tools solve this by adding a small pseudocount to unexpressed samples, but this will still result in a very high fold change. So far for the big fold-change.

If I had to guess at the reason why this is significant I would think in terms of dispersion estimation. This Bioconductor post has some discussion with the authors of DESeq2 on how this works.

Finally: calculating differential expression on these genes in these samples will never be really informative, as they really are only expressed in one sample. The way I see it you have two options:

Exclude the sample from your analysis. Straightforward, but does not solve this problem for other samples.
Filter your genes before calculating differential expression, i.e. remove all genes that are not expressed in a large number of samples. This approach is standard in for instance sleuth. A good starting point would be to remove genes that are not expressed in 75% of your samples.


Output the linear predictor from a stratified cox model? I have fitted a cox model on a pooled dataset of multiple studies, say studies A, B, and C. As these studies have different baseline hazard functions, I stratified for 'study' in the model.
Like this:
df$SurvObj <- with(df, Surv(event_rd, event == 1))

fit <- coxph(SurvObj ~ cov1 + cov2 + cov3 + cov4 + strata(study),data=df)

Now, I want to assess whether I can use the beta coefficients from the above model to predict the event probability in a new study, study D.
The baseline hazard in study D is different from those in studies A, B, and C.
When I have
predict(fit,type="lp",newdata=studyD,reference="strata")

I get this error:

Error in model.frame.default(data = studyD, formula = ~hba1c + sbp +
: factor strata(study) has new levels study=4

Why does R require that I match the strata of study D with those in study A, B, and/or C when I only want to output the linear predictor (type="lp")? As the linear predictor is independent of the baseline hazard function in this case.
I suppose it should be possible to extract the linear predictor for individuals in study D and then manually calculate the event probability using the baseline hazard of study D.
Does anybody know how to do this?
I could not find any question like this so maybe I am thinking the wrong way.

Unlike wikipedia, the linear predictor (LP) in coxph is made dependent from the baseline hazard function $\lambda$. For example, if a model has two predictors x1, x2 stratified by sex, the linear predictor would be like:
$$LP = log (  \lambda _{sex}) + \beta _{1}x _{1}+\beta _{2}x _{2}$$
Where the $log(\lambda _{sex})$ can be obtained by predict() with x1=x2=0.
We could calculate the LP without strata information (sex) by hand using the $\beta _{1} \, \beta _{2}$ provided by coef(fit) and a customised $\lambda _{unkown \, sex}$.
Below is a simulation comparing the results from predict() and by hand coef(fit).
## Toy data
df <- list(time=c(4,3,1,1,2,2,3),
          status=c(1,1,1,0,1,1,0),
          x1=c(0,2,1,1,1,0,0),
          x2=1:7,
          sex=c(0,0,0,0,1,1,1))

## Fit a stratified model
fit <- coxph(Surv(time, status) ~ x1 + x2 + strata(sex), df)
coef(fit)
#         x1          x2
# 0.79451881 -0.01917633

## Baseline linear predictor
lambda_sex1 <- predict(fit, newdata=list(x1=0, x2=0, sex=0))
# -0.746578
lambda_sex2 <- predict(fit, newdata=list(x1=0, x2=0, sex=1))
# -0.1497816

## by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=0))
# 0.009588167
predict(fit, newdata=list(x1=1, x2=2, sex=1))
# 0.6063845

## by hand
lambda_sex1 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.009588167
lambda_sex2 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.6063845

Finally, to calculate "a" linear predictor for an unknown sex:
# prevented by predict function
predict(fit, newdata=list(x1=1, x2=2, sex=2))
#Error in model.frame.default(data = list(x1 = 1, x2 = 2, sex = 2), formula = ~x1 +  :
# factor strata(sex) has new level sex=2

# by hand, assuming lambda for sex 2 is 0
0 + coef(fit)[1]*1 + coef(fit)[2]*2
# 0.7561661


sorting BAM file error using samtools I have few bam files and would like to get read counts using

samtools idxstats

[Data is aligned to hg19 transcriptome]. To use that command I need a sorted bam file. So to sort them I gave the following command.

samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam

This ended up showing:

[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
[E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes
[bam_sort_core] truncated file. Aborting.

And the tmp dir has 6 input.sorted.nnnn.bam files.

The cause of the error is that your input file is truncated or /tmp is running out of space. If you can do samtools view -H input.bam without error (reading the header also checks for the magic number at the end of the file) then it's like the that /tmp is filling up.


Tools to do VCF to MAF and MAF to VCF conversion? Normally, I would use the vcf2maf scripts to convert a VCF to a MAF (or vice versa).

This is great software, but on my system, perl scripts with dependencies are easy to break. (Here it uses VEP.)

Are there any other alternatives to this?


Fast and reliable alternatives to blast After some unexpected results (and previously reported) I heard that there are other tools for finding similar sequences besides blast that are faster and more accurate.

I only found hmmer, but I don't know if it is faster and more accurate.

Is there any review comparing such tools? A quick google scholar search didn't report anything useful.

I want to search with short DNA sequences (~400), mainly interested for taxonomic imputation so I'll be looking for similar sequences and not distant homologs. I won't try to model genes and the sequence will be unspliced as it is from direct sequencing of a a region of the 16S gene.

Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA


How to compare groups using WGS data? We have whole genome sequencing data for patients (not-cancer) (n=60) and for healthy controls (n=20). The sequencing centre has provided us with the best practice bioinformatics analyses including reads mapping (.BAM) and variant calling using GATK (.vcf) as well as annotation (annotated .vcf and .gVCF).

What should be our next steps? We are interested to see if there are any differences (global and/or specific) between the groups.

I'm not familiar with the program, but apparently Hail is setting itself up as a swiss-army chainsaw project for doing downstream analysis on variant-called datasets.

An overview of Hail can be found here:

http://blog.cloudera.com/blog/2017/05/hail-scalable-genomics-analysis-with-spark/

A tutorial on association testing can be found here:

https://hail.is/hail/tutorial.html#Association-testing


How to retrieve dataset names from python statsmodels? The Python statsmodels library grants access to many R-datasets which can be downloaded e.g. with:

import statsmodels.api as sm
df = sm.datasets.get_rdataset(dataname='iris', package='datasets').data

But for this you have to know the id or the name of the dataset.

sm.datasets.get_rdataset??
Signature: sm.datasets.get_rdataset(dataname, package='datasets', cache=False)
Source:
def get_rdataset(dataname, package="datasets", cache=False):
    """download and return R dataset

    Parameters
    ----------
    dataname : str
        The name of the dataset you want to download

There is a list of available datasets here, but this list does not contain the identifiers. Also not the detailed view.

Is there an overview of all the available datasets that contains the ids e.g. "iris"?

There is an index that is supposed to list all the datasets, but e.g. 'rice' as a name does not work.

Here is the website of the creator of rdatasets I guess with the info you are looking for:

http://vincentarelbundock.github.io/Rdatasets/datasets.html

Here is a short look at it:

Package     Item    Title   Rows    Cols    n_binary    n_character     n_factor    n_logical   n_numeric   CSV     Doc
boot    acme    Monthly Excess Returns  60  3   0   1   0   0   2   CSV     DOC
boot    aids    Delay in AIDS Reporting in England and Wales    570     6   1   0   0   0   6   CSV     DOC
boot    aircondit   Failures of Air-conditioning Equipment  12  1   0   0   0   0   1   CSV     DOC
boot    aircondit7  Failures of Air-conditioning Equipment  24  1   0   0   0   0   1   CSV     DOC
boot    amis    Car Speeding and Warning Signs  8437    4   1   0   0   0   4   CSV     DOC
boot    aml     Remission Times for Acute Myelogenous Leukaemia     23  3   2   0   0   0   3   CSV     DOC
...

I assume the Item column are the identifiers.


Retrieving Metabolic networks from WIT Database How can I get the dataset for all the 43 metabolic networks for this paper:
https://www.nature.com/articles/35036627

It looks like the WIT database does not exist anymore.

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

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

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


Getting the licenses of all conda installed packages I would like to get a list of the licenses of all packages I installed using Conda. pip has what I need in the form of pip-licenses, is there a similar thing for conda?

To my knowledge there's no single listing of the licenses, but you can:

conda list | awk '{if(NR>3) printf("%s=%s", $1, $2)}' | xargs conda info | grep license

Allegedly conda info --license should work, but at least in conda 4.6.3 it doesn't show any license information. It unfortunate that there's no convenient way to get a tabular list of package:version:channel:license information from the command line.


How to interpret fish plot in R I am working on subclonal evolution of multiple myeloma. After literature survey, I found people use Fish plot to represent tumor evolution. The example of fish plot can be found here. But I am not able to interpret this plot. How this plot represent tumor evolution and subclones in tumor. Can anybody explain with example. Thanks.

Each color is a clone (set of mutations) or an individual mutation (depends on how the input data was generated). This is a cancer with the bottleneck likely being induced by a treatment (e.g. UV, chemo, drugs, etc...)

Based on the vignette this is likely an AML (leukemia).

It looks like grey are healthy cells. At day 0 a specific clone (red) existed at a low frequency (~20-40%). The patient (or cells) was then exposed to a treatment and on day 150 another measurement was taken. In this case two novel clones (subclones of red) have taken over the population. They seem to co-exist or are co-dominate clones with orange being 50-60% of the population and yellow being 30-40%.

The example says "Simple plot with multiple subclones" so: the orange and yellow populations have additional (likely mutually exclusive) mutations that were not present in the red clone - but biologically also have all the mutations in the red clone. So my first assumption was right - each color is a proxy for a set of mutations.

You could have generated the mutation data (i.e. allele frequencies) from a wide variety of sources (scRNA-Seq, exome-sequencing, PCR genotyping, etc...)


mBED has time complexity $O(n \lg n)$, claimed by Clustal Omega paper, why? In the Clustal Omega paper, it says the following about mBED ($N$ is the number of input sequences, I assume):

we use a modified version of mBed (Blackshields et al , 2010), which has complexity of $O(N \log N)$

But in the mBED paper it says the following (section 2):

A number of $t$ sequences are initially sampled from the input dataset $X .$ Following the LLR
algorithm, this value is set by default to $t=\left(\log _{2} N\right)^{2} .$

After the seed sequences in $R$ have been chosen, all sequences in the input data are associated with a $t$-dimensional vector. This is done simply by computing the distances from all sequences to each of the $t$ seeds.

Then my conclusion is that this would lead to a complexity of $O(N (\log N)^2)$ because for each  sequence (in total $N$ of them), distances to $t = (\log N)^2$ sequences must be calculated. Where did I do wrong?

Turns out that I was dumb: I misread.
The Clustal Omega paper (in the supplementary materials) chose $t = \log N$, which I think originates from this claim from this part of the Linial et al. paper:

In random polynomial time $X$ may be embedded in $l_{2}^{O(\log n)}$ with distortion $\leq(1+\epsilon) \cdot c_{2}(X)$ (for any $\epsilon>0)$

Several questions still remain though:

The newest implementation of Clustal Omega (that I downloaded) has source code that instead chooses $t = (\log N)^2$, why changed?
The original mBED paper chooses $t = (\log N)^2$, but it was not justified on how they chose this number instead of $\log N$


Error while loading a big file with Cytoscape I installed Cytoscape 3.7.1 with java 1.8.0_191 in windows server 2012. I have a 7.4 GB csv file (about 1,500,000 reords) and when I tried to load it into the Cytoscape it throws an error java.lang.OutOfMemoryErroe:Java heap space.

My system has 500GB ram and then I changed Cytoscape.vmoptions file according to -Xms102400M but when I load my file, in task manager, the memory increased until 50GB and then Cytoscape crashed and again it threw the same error: java.lang.OutOfMemoryErroe:Java heap space.

My data is like this sample:

45677888,test1,3453453453,test2,3534523453,235412352,3453452345,test3,235423452,test4,1980/09/02,14:13,23523525,test5
45234288,test11,34234553453,test12,353434534553,23453452352,342422345,test13,23456543452,test14,1980/04/01,14:12,2323423425,test15
4243424888,test12,3235253453,test22,3533456343453,233534532352,2343452345,test33,23345342,test44,1980/03/4,14:11,23674575,test55

My file is in CSV format. col2 is source node, col1 and col3 are attributes for source node,col 4 is destination node, col5,col6,col7 are attributes for destination node, col8 until end are attributes for edge. I have about 200,000 nodes and 230000 edges.

How can I read this file in Cytoscape?

The first thing you can try is to separate the network from the meta information. Assuming you have access to a *nix machine, you can run:

awk -F"," '{print $2,$4}' file.csv > network.csv

That will produce a file like this (based on your example data):

test1 test2
test11 test12
test12 test22

That should radically decrease the file size and will give you the best chance of loading it in Cytoscape. You can then similarly extract the attributes:

awk -F"," '{print $2,$1}' file.csv > source.column1.attributes
awk -F"," '{print $2,$3}' file.csv > source.column3.attributes
awk -F"," '{print $4,$5}' file.csv > target.column5.attributes
awk -F"," '{print $4,$6}' file.csv > target.column6.attributes
awk -F"," '{print $4,$7}' file.csv > target.column7.attributes
awk -F',' '{print $2,$4,$8}' file.csv > edge.column8.attributes
awk -F',' '{print $2,$4,$9}' file.csv > edge.column9.attributes
awk -F',' '{print $2,$4,$10}' file.csv > edge.column10.attributes
awk -F',' '{print $2,$4,$11}' file.csv > edge.column11.attributes
awk -F',' '{print $2,$4,$12}' file.csv > edge.column12.attributes
awk -F',' '{print $2,$4,$13}' file.csv > edge.column13.attributes
awk -F',' '{print $2,$4,$14}' file.csv > edge.column14.attributes

Now, try loading network.csv into Cytoscape and then load each of the node and edge attributes separately.

However, a network with 200,000 nodes and 230,000 edges is just too large to be able to manipulate in any useful way in a GUI. It will just look like a horrible hairball. I really urge you to do some preprocessing and pare the network down to something more manageable. I used to work with a human PPI network of ~13000 nodes and ~15000 edges. While I could load that into Cytoscape, that was already far too big to be able to do anything useful with. I could use the data analysis aspects of Cytoscape, but not really the visualization.


how can i calculate p value of miRNA raw reads counts using R? working on lung cancer project and need to find out which genes are deferentially expressed and which are not. for that i a have downloaded miRNA data set from GEO id is GSE110907. And the file i got has to variable gene_id and read counts. from that i have to derive p value from that data/file.

I would start by re-formulating the question in this post: Your end-goal seems to the differential expression analysis (DEA) of your samples (and you have to have to levels for this like treatment-control, mutant-WT, etc ...). P-values are context-dependent and should be calculated based on your research question/hypothesis.

DEA has been addressed by many researchers and R has a lot of packages to help you with this; DESeq2, edgeR, limma, etc ...

I would go to the Biconductor website, search for one of the three above, and go through their rich vignettes and would come back here for specific questions.


Normal score transformation This question is related to a normal score transformation which is performed in the paper by
RUAN, Quansong, et al. Local similarity analysis reveals unique associations among marine bacterioplankton species and environmental factors. Bioinformatics, 2006, 22.20: 2532-2538.

They claim that given a variable X with observations $x= (x_1,\ldots,x_n)$ that are not normally distributed, normal score transformation is performed by first creating a rank vector observations $x$, $R^x = (R_1^x,\ldots, R_n^x)$ and then using and then transforming to $x'=(x_1',\ldots,x_n)$, by $x_i\ = \phi_{-1} (\frac{R_i^x}{n+1} )$ where $\phi_{-1}$ is the inverse of cummulative normal distribution.

It's intuitively clear (but maybe a bit harder to prove?) that resulting variable $X'$ has observations normally distributed, but I have a question about time series dependence. If we had that observations $x$ were correlated, say, over 4 adjacent observations, , i.e. $x_{n+1} = f(x_{n},x_{n-1},x_{n-2})$, how does this change after the normal transformation? Can we still find a function g, so that $x'_{n+1} = g(x'_{n},x'_{n-1},x'_{n-2})$ would hold? Since we're transforming based on the rank of data, I expect this to be not possible?

As you suspected, correlations will be heavily changed by such a transformation. The reason for this is also exactly as you described, namely that using ranks rather than values for the transformation eliminates all of the information needed to actually preserve these relationships. As an example, suppose we have a single observation with rank 1 across all observations that none the less shows a decrease in absolute value across them (i.e., there's a correlation with a negative slope over time). Such a transformation will result in a correlation of 0.

One might be able to find a different transformation that preserves this relationship, or might be able to simply pool all of the time points together before computing the ranks (this would then at least partially preserve some correlations, though depending on what's needed to properly normalize the data it may cause other problems), but in general there's no guarantee that this is possible.


How can I improve the yield of MinION sequencing runs? This is a frequently-asked question within the nanopore community. Oxford Nanopore currently claims that they are able to generate run yields of 10-15 gigabases (e.g. see here and here), and yet it's more common to see users only managing in the 1-5 gigabase range.

So why the big difference in yield?

I attended a talk by Josh Quick at PoreCampAU 2017, in which he discussed some common barriers to getting both good sequencing yield and long read length. It mostly boils down to being more careful with the sample preparation. Bear in mind that the MinION will still sequence a dirty sample, it will just be at a reduced yield. Here are my notes from that talk:

The hardest thing about MinION sequencing is getting the sample in the first place
There are lots of different sample types and extraction methods
You can't get longer reads than what you put in; shit in leads to shit out
DNA is very stable when not moving, but very sensitive to lateral damage
The phenol chloroform method of DNA extraction is very good, and can be used with a phase-locked gel to make extraction easier
A simple salt + alcohol extraction might be the best method for extraction (because it involves the least amount of work on the DNA)
EDTA (e.g. as found in TE buffer, and many extraction kits) is not compatible with the rapid kit
The most consistently good Nanopore runs produced by Josh's lab were 1D ligations runs on R9.4; the best overall run was a phenol-chloroform extraction + rapid kit
John Tyson can tune himself out of a low-yield hole (via software)
Getting small numbers of short reads is very important
Suggested (and mostly untested) purification techniques: spin column (60-100kb); ethanol extraction (100-150kb), dialysis (150-250kb); low melting-point agarose plug (~1Mb, DNA extraction in situ)
The nanopore protocol input is in nanograms, but should really be stated as molarity; the kit expects about 0.2 pmol input
Picture molecules tethered to the surface of the membrane. You can then see that the flow cell density is independent of the sequence length
Tapestation, Qubit and Nanodrop are all a good idea; a DNA sample that can pass all three tests will work well: no short-length shoulders by Tapestation, sufficient DNA by Qubit, high purity by Nanodrop
RNA can interfere with sequencing; digesting RNA is recommended
Freezing DNA is a really bad idea. The ice crystals are very good at chopping DNA up into small pieces
DNA that is kept in the fridge is remarkably stable; can be kept for over two years (and probably indefinitely)

Update on yield:

We had a run in June 2018 which produced about 3 million reads from PCR-amplified mouse cDNA (read N50 of ~1kb), and that was in a situation where the it was the first sample prepared by someone, the input DNA amount was probably about 20% of what it should have been, and we suffered a dataloss incident. If that's what happened with a bad run, I have high hopes for our runs in the future.

Our next cDNA run in August 2018 produced about 7 million reads and a similar read N50 (i.e. total yield of about 7Gb). The run was interrupted by a Windows update during the first night (I had disabled automatic updates previously, but they had been re-enabled), and there have since been additional software and hardware updates since to improve sequencing yield. If we have a good flow cell and similar sample prep, I expect that the next cDNA run that we'll do should produce over 8 million reads, and possibly over 10 million.


Can the canu assembler output a fastq file of the final assembly just like HGAP4? I have assembled some genome from Sequel PacBio data both with HGAP4 on the SMRT Link interface and using canu on the command line. The HGAP4 assembler outputs a fastq file of the final assembly such that I can see how reliable the sequence is along the genome. On the other hand, canu does no such thing, which means I do not know how much I can trust each nucleotide in the assembly? Is there a way canu can output the fastq file as well or at least some depth information at each base?

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

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


Is there a standard definition for "assembly polishing"? Is there a standard definition for "assembly polishing" in the field?
Is there a standard definition for what polishing algorithms do?

My understanding of "polishing" is strongly influenced by Pilon:

Whereby users input a reference FASTA along with the BAM of aligned reads against this reference FASTA. The polishing algorithm then does alignment to discover local misassemblies and other alignment inconsistencies. These are then "corrected", i.e. "polished".

In long-read assembly, "polish" refers to the step to improve the base accuracy of contig sequences. I believe the terminology was originated from the HGAP paper:

  The final consensus–calling algorithm Quiver, which takes into account all of the underlying data and the raw quality values inherent to SMRT sequencing, then polishes the assembly for final consensus accuracies in excess of 99.999% (QV of >50)...

Pilon was published more than a year later and in fact, neither its paper nor README mentions word "polish".

By the way, "error correction" refers to the step to fix sequencing errors in raw reads, not in contigs. In Canu/HGAP/Falcon, that step is done by pbdagcon or falcon_sense.