Plot to show the expression of genes between tumor and normal I have RNA-seq raw counts data for 50 samples. 20 Normal and 30 tumor. After differential analysis I got 30 DEGs. I want to make a violin plot showing the expression of each gene. I transformed counts to logCPM.

counts:

Genes   Tumor1  Tumor2  Normal1 Normal2
RP11-351J23.1   0   5   6   0
MIR17HG 989 896 8   0
RP11-563N12.2   0   0   0   0
LINC04393   0   2   16  0
RP11-336A10.4   0   0   0   0
DRH6-AB1    53  13  39  9
RP11-115J16.1   0   0   50  6
LINC70518   2   65  0   0

logCPM <- cpm(counts, prior.count=2, log = TRUE)

                Tumor1   Tumor2  Normal1  Normal2
RP11-351J23.1 11.84477 13.09301 15.71337 11.84477
MIR17HG      19.84847 19.79600 16.10350 11.84477
RP11-563N12.2 11.84477 11.84477 11.84477 11.84477
LINC04393     11.84477 12.47723 17.06532 11.84477
RP11-336A10.4 11.84477 11.84477 11.84477 11.84477
DRH6-AB1      15.72257 14.03896 18.32772 19.19283
RP11-115J16.1 11.84477 11.84477 18.68262 18.61229
LINC70518     12.44599 16.08368 11.84477 11.84477

I want to make violin plot for each gene showing expression for Tumors and Normals with this.

I first transformed logCPM. then tried like below:

ggplot2.violinplot(data=logCPM, xName='RP11-351J23.1',
    groupName='Tumor', position=position_dodge(0.8),
    backgroundColor="white", groupColors=c('#999999','#E69F00'),
    legendPosition="top")

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

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

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


Parsing gene expression value for a list of genes from GEO I am looking for the gene expression data in GEO GSE15543 with the following keywords,
(Pancreatic Beta Cells) AND "Homo sapiens"[porgn:__txid9606]

I would like to parse the gene expression value of a list of genes from all the samples reported in any particular study.I understand there are two different formats for downloading datasets, SOFT and MINiML.
I have downloaded the soft.gz file and I am using BioConductor to read the contents from soft file using the instructions given here . The file has the content that I 'm looking for , I have checked the contents on the geo2r webpage. I'm still facing challenge in finding out how to fetch out the gene expression values for a list of genes. Any suggestions?
From the files that I download, I wish to parse the gene expression value and the Sample Title and Sample Description so that I can classify the values based on disease state and control state. Any suggestions on which format can be downloaded for the above-mentioned task?

GEO provides a GUI for this (just click on the "Analyze with GEO2 button"), you can even check/download the R/Bioconductor code used for the analyses.
Basically you will create groups of samples to compare/contrast and R/Bioconductor will run behind the scenes. Here is a nice tutorial.


Techniques for analyzing and quantifying sample bleed through in genotyping with Illumina I am looking at the presence of viral genotypes within individual samples within an assay. Often times there is a sample whose read counts are firing off the charts and this sample tends to "bleed through" to the other samples.

I have recently tried confronting this problem with chi square analysis of conditional probability. For instance let's say I observe with 10 sample viral genotypes A and B. I then calculate the expected number of samples having both genotypes A and B from the number of A positive samples and number of B positive samples.

So far this is all I have. I do not know if anyone else has any ideas on how to determine and quantify bleed through in viral genotyping.

Slightly more information about this situation: these samples were run on Illumina HiSeq, we are planning on doing another HiSeq run with even more samples shortly so would like to take proper precautions with that, additionally the current pipeline utilizes Novobarcode as a demultiplexer so if anyone has suggestions on what I can implement during processing to minimize the effect of cross contamination that'd be helpful.

Are you using dual-index barcodes with different barcodes at each end? There is a known phenomena of "index switching" that occurs in Illumina reads. One way to control for this is to add in unamplified water samples that have unused barcode combinations where each of the indexes is shared with another sample.

More information can be found on this SeqAnswers thread. Here are some useful associated links:

Signal-spreading preprint -- detailed investigation of the phenomenom, including methods, graphs, and associated data
Crossblock -- tool by Brian Bushnell that attempts to remove cross contamination
Summary by James Hadfield -- mentions a few fixes suggested by Illumina, many centred around not using barcode combinations that could conflict
Illumina's quick-fire response -- they've known about it for a long time, and point out that it shouldn't be a problem for most purposes

I have suspicions that at least some of the index switching is being caused by ligation of two separate barcoded fragments during the sample prep (which definitely happens during sample prep for nanopore), but I haven't done any experimentation to work out if this is also the case for Illumina (due to money/time/goal constraints). The experimental design would be to take a prepared Illumina library just prior to sequencing, ligate on ONT adapters, then sequence on both the ONT and Illumina machines to see if chimeric reads discovered via nanopore sufficiently matches the index switching observed via short-read sequencing.


parsimony and maximum likelihood tree comparison in R Maximum Likelihood analysis compares the tips or the species based on their sequence similarity while Parsimony analysis compares the characteristic features among the species.

So far I did comparison of both the tree using dendextend. Is there any other methods or library to make the comparison more quantitative.Due to my limited understanding I have this in mind. Can i even compare these two methods?

Can i even compare these two methods?

See if this paper helps:
https://www.nature.com/articles/nature02917


knnImputation for missing data I have a dataframe containing missing values. I want to impute them with knnImputation as shown below.

Now my question is: are data imputed row wise or column wise and is there a way to give this option in knnImputation?

A <-c(1,2,3,4,5, NA, NA)
B <-c(2,3,NA, 2,4, 7, 6)
C <-c(2,3,4, NA, 5, 1,7)
D<-21:27
data<-data.frame(A,B, C, D)
new_data<-knnImputation(data, k = 3, scale = T, meth = "weighAvg", distData = NULL)

The data are imputed using values from the same column

Here's a simplified example to show what's going on:

A <-c(1,2,3,4,5)
B <-c(2,3,NA, 2,4)
C <- c(2,3,4, 6, 1)
data<-data.frame(A,B,C)
data

#>  A  B C
#>1 1  2 2
#>2 2  3 3
#>3 3 NA 4
#>4 4  2 6
#>5 5  4 1

# I'm using scale = F to make things simpler
knnImputation(data, k = 3, scale = F, meth = "weighAvg", distData = NULL)

#>  A        B C
#>1 1 2.000000 2
#>2 2 3.000000 3
#>3 3 2.594272 4
#>4 4 2.000000 6
#>5 5 4.000000 1

# calculate euclidean distances, dropping column of interest
dist(data[,-2])

#>1        2        3        4
#>2 1.414214
#>3 2.828427 1.414214
#>4 5.000000 3.605551 2.236068
#>5 4.123106 3.605551 3.605551 5.099020

# use k rows with smallest distance to row of interest to calculate weighted average
# using exp(-dist) as weights
weight_sum <- weight_sum = exp(-2.828427) + exp(-1.414214) + exp(-2.236068)
2*exp(-2.828427)/weight_sum + 3*exp(-1.414214)/weight_sum + 2*exp(-2.236068)/weight_sum

#>2.594271

If you want to impute using values from the same row, all you have to do is transpose the matrix:

t(knnImputation(t(data), k = 2, scale = F, meth = "weighAvg", distData = NULL))

#>     A        B C
#>[1,] 1 2.000000 2
#>[2,] 2 3.000000 3
#>[3,] 3 3.086729 4
#>[4,] 4 2.000000 6
#>[5,] 5 4.000000 1

(I had to lower k because there are only 2 extra columns)


Reading a txt in R with weird header I am trying to read a txt file in R but always I am getting

Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  line 113 did not have 10 elements

I noticed that in the header of my file I have something like

 0000775 €   g¿´�€   4_®™00000016125 13404727020 017402  0                                                                                ustar   devons01                        stlab-icgc

Even by removing these characters I am not able to read this file

This is link of my file

https://www.dropbox.com/s/rpc5jgti5aa433b/txt.txt?dl=0

Could somebody please heap me in reading this file in R?

We can skip the first row:

read.delim("txt.txt", skip = 1, header = FALSE)

Or we can get the "uncorrupted" version of the table, it looks like it is from this paper, Table S4.


Statistical approach to compare the SNP genotyping data among set of individuals So, I have the genotyping data of about 650,000 SNPs for 96 individuals. I already know the Y DNA haplogroup of these individuals, so to some extent, I have a gross understanding of their ancestry.

What would be the best way to go about doing this? From what I have been reading, maybe a PCA would be a good idea? Suppose individuals 1-10 belong to a particular population (unique ancestary) then their SNP genotypes should be significantly similar to each other than say individuals 11-20 who belong to another population group. So, in theory I should get two clusters in my PCA plot.

I have the genotype of each SNP for each individual. How would I go about creating a PCA plot or any other alternative way of analysis for this kind of questions.

PCA sounds like a good start to see a rough population structure among your 96 individuals, this step is half a sane thing to do, but also it can give you a hint how many distinct populations there are in your dataset. I have not personally run PCA on SNP data, but it seems that there are quite lot of tutorials online. For instance this one looks  nice.

Once you know how many categories you can expect in your population, you can use Structure to classify individuals to the categories.


Why is the UniProt REST API returning multiple results, when I am only providing one ID to be checked? Can some explain why when I run a REST API query on UniProt for one ID ('Q9UJL9'), the results for three IDs are returned?
The query I'm using is:
https://www.uniprot.org/uniprot/?query=Q9UJL9&sort=score&columns=id,go(molecular%20function)&format=tab
(I didn't hyperlink it since you need to see the query to see where I'm going wrong).
When I run the query above, the output is:
Entry   Gene ontology (molecular function)
Q9UJL9  DNA-binding transcription factor activity, RNA polymerase II-specific [GO:0000981]; DNA-binding transcription repressor activity, RNA polymerase II-specific [GO:0001227]; metal ion binding [GO:0046872]; RNA polymerase II transcription regulatory region sequence-specific DNA binding [GO:0000977]
Q8WWY3  identical protein binding [GO:0042802]; ribonucleoprotein complex binding [GO:0043021]; RNA binding [GO:0003723]; snRNP binding [GO:0070990]; U4atac snRNA binding [GO:0030622]; U4 snRNA binding [GO:0030621]
Q6A162  structural molecule activity [GO:0005198]

I checked the uniprot page of the query I want (Q9UJL9), the other items are interactors (along with many others) but not synonyms/linked (it seems, e.g my query is a zinc finger, one of the linked records, Q6A162, is a keratin protein).

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


Difference between clinical significance and statistical significance I am working on a clinical data analysis project and am from IT domain.

While I am analyzing the data, I came to know that though some of the variables are not statistically significant, they are clinically significant.

I also read that difference in median values between case and control group for variables is used as a measure to indicate clinical significance

Can someone help me understand this with a simple example as to why difference in median between groups is considered for clinical significance?

This is non-parametric statistics, the mean requires confirmation to the normal distribution (but there are exceptions), or perhaps better put it requires a fixed relationship between mean and standard deviation. Non-parametric stats always use the median and are distribution free. Wilcoxon test/ Mann-Whitney U test  are used to replace t-test. The Kruscal-Wallace test replaces the 1-way ANOVA
Finally the Friedman test replaces 2-way ANOVA but has its critiques.  The t- and ANOVA tests require normality, F-max is one approach to determine normality.


Power law distribution alpha values I did a Pearson co-expression analysis for generating networks for my tissue-specific (chondrocytes) RNA-seq data. I used R package poweRlaw to check for the power-law distribution. We got the following values:
 Xmin = 14
 goodness of fit = 0.035
 α = 4.76

My question is what does alpha (α) of more than 4 indicates?
Is it is okay to have such values?

Thank you

It means the distribution is heavily skewed and has a 'long-tail', thus the variance is much greater than the mean. In other words most values are around zero, with a small number of values having extreme values.
[
You need to check the goodness of fit because that is saying the deviation from a power-laws distribution is just significant. So it ain't quite power-laws, but the general gist of a 'long-tailed distribution' is correct.
Its a good result and a nice analysis, but you are not quite there at present (needs to be >0.05). Check the negative binomial distribution.

I've seen you data, it looks very good and the power series is not a bad explanation of the data. There are two possibilities for the goodness of fit being 0.035

The goodness of fit is low because of 'outliers' at the beginning of the distribution
An alternative distribution, such as a negative binomial distribution fits better
Check it is 0.035 not 0.35 because that changes everything.

The best way to do this is to plot the power series expected distribution onto the graph and look where the fit 'fails'.
If it is not fitting at the beginning of the distribution thats is about outliers and easily solvable. If it is not really capturing the shape then its probably a slightly different distribution.
If its point 1 then just remove the very first point in the graph and try again and check the p-value. Or take the mean of the clump of three points next to the first point and try again (or just remove those three point). Don't do any more than that and the technical word is "outlier" analysis. In your write-up you must explain the issue around any 'outliers' and that is perfectly acceptable. If it is a power-series distrupted by an outlier(s) there is formal analysis to demonstrate an 'outlier' using outlier analysis if any one ever asks you.
If its point 2. then use e.g. negative bionomial or the gamma distribution (alpha <1) and check the goodness of fit.
You want a goodness of fit > 0.05 and preferably >0.1 to be sure its the right distribution. There are formal statistics to do this, but that is the main one.
My guess is there are outliers at the beginning because your sample size is not large enough to identify the variation correctly, hence the fit is low. In 'long-tails' its easy to have a 'low sample size' because of the way they work. Thus in this case mapping a distribution is a perfect way to compensate for the sampling.


Is there a database of disease prevalence? I want to collect data on how common different diseases are in different populations (or, at least, globally). Ideally, it would give me a way of querying the database with a disease name, and would return the number of cases per N population.

I have found FIDD, but it seems to be no longer maintained. Are there any actively maintained, similar tools that have the information I am after?

For rare diseases, you can use Orphanet. You can download an xml file with prevalence from orphadata.

There is an OWL version of Orphanet called ORDO, you can browse it on the EBI Ontology Lookup Service. See for example Arrhythmogenic right ventricular cardiomyopathy and point prevalence relationships in the lower right box. There are a variety of programmatic ways to access this, including the OLS web API, SPARQL queries, or downloading the OWL and using an ontology API.


Finding proteins in a sequence I have a file with protein sequences, all on a single line, like this:

>Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*

I want to print all sequences between an M residue and a * that are at least 6 residues long, and I would like the output like this:

>contig_1.1.1
M.......*
>contig_1.1.2
M........*
>contig_2.1.1
M............*
>contig_2.1.2
M......*

Where the first number after the underscore is the number found in as the first character of the sequence (it is the reading frame) and the second number is the sequence number. and the last number is the fragment number. So, since the 1st sequence starts with;

>Contig_1 1FFPMPQ

The name will be >Contig_1.1 and each fragment found will be >Contig_1.1.1, >Contig_1.1.2, ... , >Contig_1.1.N.

After implementing the suggestion in the answer below and escaping the * in my regular expression, my current code is:

import sys
import re

myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD>Contig_2 1FFPCHNSVVYFLLLICLLKPLWSQGVVRHWYRTRCDAWLAVML*YLHVSFIPSIHFLLGRLCTLLYLDVFELYACAAGFMPCLLR*FSLGHSLFLLPLKFPATVKFFSLLGQIMNNISLMNDLQTGEKLQETH*EVEDKSALLVSIIESRCCVTLQPQDFSN*HLGKL*KSVGQKLK*VNDQITY*"
# parse record
my_protein = []
for line in myprotein:
    if not line.startswith('>'):
        my_protein.append(line)
    else:
        header = line

# build sequence
my_pro = ''.join(my_protein)

# cleave sequence
peptides = re.findall(r'(M[A-Z]{6,}\*)', my_pro)

# build a formatted list
records = ["%sNo.%d\n%s\n" % (header, idx + 1, fragment) for idx, fragment in     enumerate(peptides)]

# print list
for record in records:
    sys.stdout.write("%s" % (record))

And it prints:

 >No.1
 MVTGSCQTLVQDKV*
 >No.2
 MIPPCVFHPIYPLFTWSPLHPIISGCF*
 >No.3
 MNNISLMNDLQTGEKLQETH*

How can I change this to get my desired output format?

The * character has a special meaning in regular expressions, it means "0 or more". So when you write {6,}*, the regex engine tries to apply the "0 or more" to the {6,*} and gets confused:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3.6/re.py", line 222, in findall
    return _compile(pattern, flags).findall(string)
  File "/usr/lib/python3.6/re.py", line 301, in _compile
    p = sre_compile.compile(pattern, flags)
  File "/usr/lib/python3.6/sre_compile.py", line 562, in compile
    p = sre_parse.parse(p, flags)
  File "/usr/lib/python3.6/sre_parse.py", line 855, in parse
    p = _parse_sub(source, pattern, flags & SRE_FLAG_VERBOSE, 0)
  File "/usr/lib/python3.6/sre_parse.py", line 416, in _parse_sub
    not nested and not items))
  File "/usr/lib/python3.6/sre_parse.py", line 619, in _parse
    source.tell() - here + len(this))
sre_constants.error: multiple repeat at position 10

In order to use special characters literally, as simple strings, you need to escape them by adding a \ before them. So try:

import re
myprotein = ">Contig_1 1FFPMPQ*CCLFSSSYLPA*AIMVTGSCQTLVQDKV*CLASCHVMIPPCVFHPIYPLFTWSPLHPIISGCF*AVCMCSWIHALFVEIILSWSLSLPFTFEISSD"
protein = re.findall('M[A-Z]{6,}\*',myprotein)
print(protein)
['MVTGSCQTLVQDKV*', 'MIPPCVFHPIYPLFTWSPLHPIISGCF*']

To get the output you asked for while reading from a file, you can do:

import sys
import re

proteinFile = sys.argv[1]

with open(proteinFile, 'r') as f:
    for line in f:
        ## Split the long input line into the list 'sequences'
        sequences = line.rstrip('\n').split('>')

        for i, k in list(enumerate(sequences)):
            ## Skip the first element which will be empty.
            ## there must be a pythonic way of doing this, but I'm
            ## a Perl hacker, so I don't know it.
            if len(k) == 0:
                continue

            seqName, seq = k.split(' ')
            ## Get the frame
            frame = seq[0]
            ## Remove the frame from the sequence
            seq = seq[1:]
            # cleave sequence
            peptides = re.findall(r'(M[A-Z]{6,}\*)', seq)

            for idx, peptide in enumerate(peptides):
                print(">%s.%s.%d\n%s" % (seqName, frame, idx+1, peptide) )


PDB residue numbering The residue numbering in PDB (in my understanding) coincides with UniProtKB align. However, this might not be consistent with that in PDB fasta file. To complicate it further, UniProtKB align can be a few mutations away from the PDB sequence. Disordered/unidentified residues in structures also create gaps in numbering.
Is there a way/suggestion to create a matched residue numbering between sequence and structure? For example
12345678 (residue numbering)
ATPWQMSQ (sequence)
---W-MSQ (PDB residues) "-" for missing residue in structure

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

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

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

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

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

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

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

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

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


Question about public availability of human SNP dataset with country of origin Is there a publicly available data set for humans SNPs, preferably together with the country of origin? If so, could someone please point me towards it.
I've looked through other questions on the site and have been unable to find a similar one. I don't know how sensitive human SNP data is, so I'm not sure if this would be public.

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


Normalize ATAC-seq/Dnase-seq sequencing reads coverage signals over estimated background I'm trying to normalize the coverage signals of ATAC-seq reads against its own background using normal distribution, described in this paper

It says:

  Finally, all open chromatin coverage measurements were normalized by
  standardization to the mean and standard deviation of coverage over a
  set of 25,000 randomly selected background regions. To select
  background regions, the set of peak open regions were widened to
  20,000 bp, reduced, and subtracted from the genome assembly.
  Thereafter, 25,000 random positions were selected and widened to
  reflect the distribution of widths in the set of open peaks. Coverage
  within these background regions was then calculated, and regions with
  zero coverage were discarded (~5%). The distribution of counts within
  background regions approximated a log-normal distribution. Mean and
  standard deviation of these background regions was calculated and used
  to transform the coverage measurements for the entire genome.

Does anyone know any existing codes doing that? I know Anshul's pipeline does that but can't find the specific script for that step.

Supplementary File 5 and 6 contain the code you're looking for.


Download complete Phylogenetic Tree as JSON I would like to integrate the phylogenetic tree into a website. Does anyone know where I can download the complete tree? Best would be a JSON, which looks something like this:
[{
  "level": "domain",
  "name": "eukaryota",
  "commonNames": {
    "en": ["eukaryotes"],
    "de": ["Eukaryoten", "Eukaryonten"]
  },
  "children": [
    {
      "level": "kingdom",
      "name": "plantae",
      "commonNames": {
        "en": ["plants"],
        "de": ["Pflanzen"]
      },
      "children": [...]
    },
    {
      "level": "kingdom",
      "name": "animalia",
      "commonNames": {
        "en": ["animals"],
        "de": ["Tiere"]
      },
      "children": [...]
    },
    ...
  ]
},
{
  "level": "domain",
  "name": "bacteria",
  "commonNames": {
    "en": ["bacteria"],
    "de": ["Bakterien"]
  },
  "children": [...]
}]

Other formats and structures are ok, too. The important thing is that the level of the nodes can be retrieved. In the end I only want to have the main levels Domain, Kingdom, Phylum, Class, Order, Family, Genus and Species in the tree. (https://en.wikipedia.org/wiki/Taxonomy_(biology)#/media/File:Biological_classification_L_Pengo_vflip.svg) So if there are more nodes from other levels like Clade, Tribe or Subfamily, it must be possible to identify and remove them.
Thanks a lot!

NCBI maintains an SQL dump of their taxonomic database
You can see what that contains here
For smaller requests (depending on what you need), you might also consider the ENSEMBL REST API - although I don't know of a way to retrieve the complete tree (shoot them an email): https://rest.ensembl.org/taxonomy/classification/53399?content-type=application/json


Finding differentially expressed genes between two samples Hi, I was wondering whether we can find the differentially expressed genes between the Double-KO and the Shox2-KO. For the FindMarker() function, only see the differential expressed gene between different identities. I don't know whether we can find the differential expressed gene between only these two samples.

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

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


Calculating bit score: How do you find lambda and K? To calculate bitscore from score you can use this equation: $S' = (lambda*S - ln(K)) / ln(2)$

If I am trying to manually calculate the bitscore of an HSP of a pairwise blastn alignment, and I know the alignment score, how do I calculate the parameters K and lambda?

Is there a more simple equation for pairwise nucleotide alignments?

NCBI BLAST uses tables of pre-computed values for lambda/K for gapped protein and nucleotide comparison and their values depend on the amino acid substitution matrix or the substitution scores + gap penalties.

You can find these tables (from line 589 for blastn) in the file c++/src/algo/blast/core/blast_stat.c from the blast source package.


How to convert a Pileup file to VCF format with Hg19 alignment I would like to know how to convert a pileup file for an adna y chromosome to vcf with Hg19 alignment.

Short answer

You'll have better luck searching for the term variant calling than with the term converting. A Google search for "variant calling" gives several beginner-level tutorials on the topic. Commonly used tools for variant calling include samtools, GATK, and FreeBayes.

Long answer

Converting data files from one format to another is a common task in bioinformatics. But that is not the task you want to perform here—the term convert is misleading.

A pileup file contains reads aligned to a reference sequence. A VCF file contains a description of putative or validated genetic variant. The process of creating a VCF from alignments is commonly called variant discovery, variant prediction, or variant calling. This is a broad and active area of research interest in bioinformatics, with many available algorithms and tools, each with a different set of parameters and assumptions.

This forum is not well designed to give you a good introduction to this topic. We are much better equipped to answer specific questions you have about how a particular tool or algorithm works, or why you are observing a particular unexpected result, and so on...assuming you have put in some time to research the question yourself.


Converting gene names from one public database format to another This is a question from /u/apivan19 on reddit. The original post can be found here.

I have some proteomics data that was given to me with the UniProt gene identifiers in column 1. I've been trying to convert these to normal gene symbols using various programs, but it is proving to be difficult.

The Uniprot website does it fairly decently, but it is not able to convert all of them and then adds some unknown genes into my list.

For example, I will give it 5439 genes in UniProt notation, and will say "5420 of 5439 UniProt identifiers have been converted to 5450 gene symbols"... which is ridiculous.

I tried using David to change the symbols, but it returns them to me in some ridiculous, random order and there's no way I can sort... actually there might be but it'll take a second.

What are some of the easiest ways to do this? It's already very time consuming and am looking for simpler solutions

I tend to use Ensembl Biomart for such queries since there are APIs for various programming languages, e.g. biomaRt, and, maybe more interestingly, via a REST API (although it’s a pretty terrible one).

To translate identifiers from different databases, proceed as follows:

Choose database “Ensembl genes”
Choose dataset your desired oganism
Go on “Filters” › “Gene:” › “Input external reference ID list”

Select the chosen source database
Provide a list of IDs, delimited by newline

Go to “Attributes” › “Gene:” › untick “Transcript stable ID”

If Ensembl IDs are desired, leave “Gene stable ID” ticked …
Otherwise untick it; go to “External:”, tick your desired identifier format

Click “Results” at the top left. This gives a preview that can be exported into various formats; alternatively the top centre buttons “XML” and “Perl” provide the query in XML (for SOAP/REST requests) and as a (horrendously formatted) executable Perl script.


combining the use of workdir and option --jobscript in snakemake It seems that in snakemake the script specified after the --jobscript cannot be used properly if a workdir: is specified in the snakefile.
The path of the script specified becomes relative to the workdir defined in the snakefile instead of being relative to the current working directory. I do not know if this is a feature or a bug !!

Nevertheless, it becomes quite painful to workaround this problem. The only way I found was to parse the command line in the snakefile, for example:

snakemake -c "qsub" -j 30 --js ./sge.sh --latency-wait 30 -rp

in order to copy ./sge.sh in the output directory defined by the keyword workdir in the snakefile.

As a side effect, if you specify option for qsub on the command line, for example:

snakemake -c "qsub -e ./logs/ -o ./logs/" -j 30 --js ./sge.sh --latency-wait 30 -rp

the logs folder must be created under the workdir directory.

Is there something I don't understand with the --jobscript option? Am I not using it the right way? Am I going against the best practice?

To quote from the snakemake documentation

All paths in the snakefile are interpreted relative to the directory snakemake is executed in. This behaviour can be overridden by specifying a workdir in the snakefile:
workdir: "path/to/workdir"
Usually, it is preferred to only set the working directory via the command line, because above directive limits the portability of Snakemake workflows.

It's not explicitly stated, but it's sort of implied that all relative paths then become relative to the working directory. I would expect that specifying an absolute path would get around this.
As an aside, in my mind setting the working directory is usually only needed on clusters without a shared file system (presumably shared between the worker nodes, but not with the head node), since there you can't cd $WORKDIR before running snakemake. This is then normally done with your scheduler, in such cases.


How to find common sequences among 6 multi-fasta files I have 6 multi-fasta files, every of them contains ca 1500 sequences like that:

>Haladaptatus sp. R4
MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
ALEEVADR
>Methanosphaera sp. WGK6
MAKVNVYSLKGDITEEIELPEIFEEEYRPDVIKRAVISTQTARIQPWGANPMAGKRTTAE
SFGSGRGAAMVPRVKGSSKAAFVPQAIGGRKAHPPRVNTIYHEKINKKERILAIRSAIAA
TANKEIVEQRGHAVANLEQVPFVVDDELETIKTTKETREIFKDLGIMDDILRAKKGRKIK
SGKGKLRGRKYRTPKGPLVVVGNDRGISLGARNHAGVEVVEVNNINAELLAPGTHAGRLT
IYTKSAVEKLADLFQQNRS

and so on.
I need to create a new result file containing sequences that are common in these multi-fasta files.
Example:
If sequence

>Haladaptatus sp. R4
    MQATVRDLNGEDADTVDLPDVFETTVRTDLIKRAVLAAQANRKQDYGTDPHAGMRTSAES
    PGSGRGMAHVPQTNGRGARVPFTVGGRVAHPPKAEKDRSRSINKKERKLAVRSAIAATTD
    AERVSERGHRFDEDTELPLVVSDDFEDLVKTQEVVSFLEAVGIDADIARAEDNKKVRAGR
    GTTRGRKYKTPKSILFVTSEEPSRAARNLAGADVATAREVNTEDLAPGTQAGRLTVWTES
    ALEEVADR

occures in 1 and 2 and 3 and 4 and 5 and 6 file it shoud be written to the result file.
May anoyone help?

Edited: this is my Python script which should work on 2 files (not 6 for now to make it simplier):

for z in f1:
if z[0] == '>':
    znaleziony = False
    gatunek1 = z[1:len(z)]
    #print(gatunek1)
    for x in g1:
        #print(x)
        if x[0] == '>':

            gatunek2 = x[1:len(x)]
            if gatunek1 == gatunek2:
                wynik.write('>')
                wynik.write(gatunek1)
                znaleziony = True
                continue

elif znaleziony == True:
    wynik.write(z)

wynik.close()

Unfortunately it repeats some identifiers of sequence and does compare only identifiers (not whole sequence).

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

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

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

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

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

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

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

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

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

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

Finally, convert to fasta again:

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

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

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

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

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


CNSeg no method coercing S4 class to vector I'm attempting to analyze some genomic data with the help of a package called CNSeg. The package is located here and the documentation is here.

My input data looks like:

> str(segments)
'data.frame':   11897 obs. of  7 variables:
 $ X                : int  0 1 2 3 4 5 6 7 8 9 ...
 $ SegmentID        : int  72 73 74 75 76 77 78 79 80 81 ...
 $ Chromosome       : int  1 1 2 2 2 3 3 3 3 3 ...
 $ StartPosition    : int  754192 145260908 21494 141215321 141275624 63411 69812903 69884262 126473310 126790130 ...
 $ StopPosition     : int  145258178 249212878 141214996 141275051 243052331 69811900 69884106 126457276 126772699 197852564 ...
 $ Median.Log2.Ratio: num  -0.014 0.311 -0.003 0.059 -0.012 -0.018 -0.106 0.007 -0.171 0.001 ...
 $ FileName         : Factor w/ 95 levels "TSB02183","TSB02184",..: 1 1 1 1 1 1 1 1 1 1 ..

(Note: I did import the segments data with FileName as a vector instead of a factor as well. No difference in the error below).

I followed the vignette to the letter and, at the end of the vignette, the creator tells me that I can use the dist() method. When I try to use it, I get an error: Error in as.vector(data) :
  no method for coercing this S4 class to a vector.

Here is my session information:

> segments <- read.csv("Probe_Segments_CN.csv")
> cnseg <- CNSeg(segList = segments, chromosome = "Chromosome", end = "StopPosition", start = "StartPosition", segMean = "Median.Log2.Ratio", id = "FileName")
> rdseg <- getRS(cnseg, by = "region", imput = FALSE, XY = FALSE, what = "mean")
Processing samples ... Done
> data("geneInfo")
> geneInfo <- geneInfo[sample(1:nrow(geneInfo), 2000), ]
>
> rdByGene <- getRS(cnseg, by = "gene", imput = FALSE, XY = FALSE, geneMap = geneInfo, what = "median")
>
> reducedseg <- rs(rdseg)
> f1 <- kOverA(5, 1)
>
> ffun <- filterfun(f1)
>
> filteredrs <- genefilter(rdseg, ffun)
> filteredrs <- madFilter(rdseg, 0.8)
> dist(filteredrs)
Error in as.vector(data) :
  no method for coercing this S4 class to a vector
> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] tools     stats     graphics  grDevices
[5] utils     datasets  methods   base

other attached packages:
[1] CNTools_1.34.0    genefilter_1.60.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17         AnnotationDbi_1.40.0
 [3] BiocGenerics_0.24.0  splines_3.4.0
 [5] IRanges_2.12.0       bit_1.1-14
 [7] lattice_0.20-35      xtable_1.8-2
 [9] blob_1.1.1           parallel_3.4.0
[11] grid_3.4.0           Biobase_2.38.0
[13] DBI_1.0.0            survival_2.41-3
[15] bit64_0.9-7          digest_0.6.15
[17] Matrix_1.2-9         S4Vectors_0.16.0
[19] bitops_1.0-6         RCurl_1.95-4.10
[21] memoise_1.1.0        RSQLite_2.1.1
[23] compiler_3.4.0       stats4_3.4.0
[25] XML_3.98-1.11        annotate_1.56.2

If you are using the latest version of Bioconductor you should always look at the latest documentation of the package in Bioconductor. There could be changes from 2011 to 2018 (In this case there aren't any changes in the section you mention).

The problem is the package doesn't define how dist should behave with that object class (filteredrs is a RS object if I read the source code well).

You can extract the metrics to calculate the distance yourself or ask the maintainer to add this method to the class. However it could be that the maintainer answer you if you post the question in the support site of Bioconductor.


What does it mean when PDB lists multiple organisms for the same structure? I'm browsing through ribosomes which have been modeled with Cryoem on PDB and am pretty confused by the fact that for some, PDB, or rather the authors of the deposition probably, put multiple organisms from different domains of life... How does one interpret that?

So far as i understand bos taurus is a cow and e.coli is a bacterium. So eukarya and bacteria.. Thank you.

Since this is a 40S ribosomal protein, and therefore eukaryotic-specific (as the description says) the listing of E. coli is not intuitive as you note. If you go down to "Macromolecules" and click on "Nucleic Acids/Hybrids" you can see the listed E. coli sequence is just 28bp, and all Us (poly(U) tail?). So, you can safely ignore the bacterial annotation I think.


Normalization for two bulk RNA-Seq samples to enable reliable fold-change estimation between genes I have two bulk RNA-Seq samples, already tpm-normalized.

I would like to know what is a reasonable normalization procedure to enable downstream log fold-change estimation.

The distribution of the two samples using the common set of genes looks similar:

However, the two samples have only been tpm-normalized, is it enough to guarantee reliable fold-change estimation? Should I use another normalization procedure, e.g. Quantile Normalization, before comparison?

My objective is to define a signature using the genes that are up-regulated in Sample1 with respect to Sample0, and vice versa. I'm using log fold-changes, but I'm concerned that their value may be affected by each sample distribution.

Do you also have suggestions for the definition of up-regulated genes with these data?

You have only two samples?

You aren't going to be able to draw strong conclusions from that no matter what you do.  Clever statistics don't work without replicates.


Hba-a2 has two copies on UCSC genome browser I was browsering the UCSC viewer of the genome at this region.

UCSC genes label two genes as Hba-a2 while GENCODE thinks the left one is Hba-a1

Which annotation is correct?

Gencode is correct (when in doubt, always assume Gencode/Ensembl is correct when it differs from UCSC). What presumably happened is that the UCSC annotation pipeline noticed that Hba-a1 and Hba-a2 are nearly identical (they produce the same protein and have minimal nucleotide difference) and just misannotated.


How to create a graph of multi column gene dataset? well, I'm new to Bio. what I'm going to do is create a network/graph of a dataset related to the genes I think. here's the sample of the nodes dataset with CSV format named lung-nodes :

and this is a sample of edges I think with trrust name :

these data are stored in two separate file

what I must do is create a graph(network) of these nodes and edges by Python and a visualization library like networkX, but I do not how.
I must make a relation between these to the dataset by creating a graph of them
can you help me out?

As far as I understand, it seems that the first table you provided contains biological/experimental information related to your genes, i.e. the nodes of the network you will build. The second table is the one you need to build the network: each line is an edge. You can use it as input and build the network using the IGraph library which works for both Python (as you requested) or R (as you suggested as an alternative).

Most probably, but I can't say for sure, the first table was just a dataset from which you (or someone for you) recovered a list of nodes that where then submitted to a database, i.e. String, commonly used to recover information related to the interactions between your nodes, i.e. genes. Once you have this info (in the second table) you can build the network!

edit:

if you're not familiar with IGraph, you can have a look at this question describing how to obtain an IGraph graph from a csv file or this other one which uses a different csv file format to build a graph


Filtering out all seqs with mutations of list2 from list1 I have 2 lists, list1 and list2, of protein sequences of the same given gene in different strains. In list1 are reported amino-acid sequences from resistant strains, and in list2 are reported amino-acid sequences from susceptible strains. Each item of each list is an amino-acid sequence.

I want to filter out (in Python 3) from list1 all of the variants that are also present in list2. The final aim is to have a list of the sequences that have unique mutations that are not in the 2nd list. The mutations can be of any type (e.g., indels, mismatches).

All the sequences are of at least 90% identity to some reference sequence  (filtered out the low matches using seq alignment).

Any enlightenment / package / tool for this job will be appreciated!

Toy example:

a, b, c : amino-acids

For this input:

Reference seq: aaabbb
Resistant strains: list1 = [aacbbc, ababbb]
Susceptible strains: list2 = [aacbbb, aaabbc]

The expected output would be: [ababbb], since the mutations in list2 cancels out the item aacbbc.

Of course this is a simple example and it is possible to have more complex scenarios (e.g., longer indels and mismatches).

Samtools/VCF is overkill for comparing amino acid sequences from a single protein. I would suggest making pairwise alignments between each sequence in list2 and the reference and making a list of mutations (2c and 6c in your example). Then repeat the procedure for list1, discarding sequences if they contain mutations in your stored list. This can all be done using biopython. Here's how to make the alignments to get you started:

from Bio import pairwise2
from Bio import SeqIO
seq1 = SeqIO.read("alpha.faa", "fasta")
seq2 = SeqIO.read("beta.faa", "fasta")
alignments = pairwise2.align.globalxx(seq1.seq, seq2.seq)


Random addition method for phylogenetic tree reconstruction I have been working on nucleotide and protein data to reconstruct maximum parsimony trees using MEGA software to get parsimony tree, via "random addition" method is used to generate the initial tree(s).

I need to know that when randomly chosen sequences are added to reconstruct an initial tree, on what basis the best tree is chosen??

I would appreciate if anyone can guide me about this method or suggest me some good literature about this in order to make me understand the concept of this method.

MEGA manual may help to explain the criteria used in selecting the final maximum parsimony tree.

The random tree additions refers to the generation of random initial trees for a heuristic search to find MP trees. In this case, a tree is generated by randomly selecting a sequence and adding it to the growing tree on a randomly-selected branch.

Instead of using the random tree additions you can use the min-mini algorithm for the initial tree set.


How to do `bedtools intersection` using pandas alone? I have two pandas Dataframes, using python3.x:

import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)
print(df2)
##       0    1    2         3  4  5
## 0  chr1  155  200  feature5  0  -
## 1  chr1  800  901  feature6  0  +

The columns to focus on in these dataframes are the first three columns: location, start, and end. Each start:end value represents a distance on location (e.g. chr1, chr2, chr3).

I would like to output the intersection of df1 against df2. Here is the correct output:

chr1    155 200 feature2    0   +
chr1    155 200 feature3    0   -
chr1    900 901 feature4    0   +

What is the most efficient (in terms of runtime and RAM) to find the intersections using pandas?
There are python wrappers around BEDTools: https://daler.github.io/pybedtools/intersections.html

I was hoping there was a pandas only solution...

Your setup:

import pandas as pd

dict1 = {0:['chr1','chr1','chr1','chr1','chr2'],
    1:[1, 100, 150, 900, 1], 2:[100, 200, 500, 950, 100],
    3:['feature1', 'feature2', 'feature3', 'feature4', 'feature4'],
    4:[0, 0, 0, 0, 0], 5:['+','+','-','+','+']}

df1 = pd.DataFrame(dict1)

print(df1)

##       0    1    2         3  4  5
## 0  chr1    1  100  feature1  0  +
## 1  chr1  100  200  feature2  0  +
## 2  chr1  150  500  feature3  0  -
## 3  chr1  900  950  feature4  0  +
## 4  chr2    1  100  feature4  0  +

dict2 = {0:['chr1','chr1'], 1:[155, 800], 2:[200, 901],
    3:['feature5', 'feature6'], 4:[0, 0], 5:['-','+']}

df2 = pd.DataFrame(dict2)

Now use pyranges (https://github.com/biocore-ntnu/pyranges):

# pip install pyranges
# conda install -c bioconda pyranges

import pyranges as pr
cols = "Chromosome Start End Name Score Strand".split()
# Add column names to the dfs so PyRanges knows which
# represent Chromosomes, Starts, Ends and optionally Strand

df1.columns = cols
df2.columns = cols

# create PyRanges-objects from the dfs
gr1, gr2 = pr.PyRanges(df1), pr.PyRanges(df2)

# intersect the two
gr = gr1.intersect(gr2)

# get the result as a   df
df = gr.df

Results:

print(gr)
# +--------------+-----------+-----------+------------+-----------+--------------+
# | Chromosome   |     Start |       End | Name       |     Score | Strand       |
# | (category)   |   (int32) |   (int32) | (object)   |   (int64) | (category)   |
# |--------------+-----------+-----------+------------+-----------+--------------|
# | chr1         |       155 |       200 | feature2   |         0 | +            |
# | chr1         |       900 |       901 | feature4   |         0 | +            |
# | chr1         |       155 |       200 | feature3   |         0 | -            |
# +--------------+-----------+-----------+------------+-----------+--------------+
# Stranded PyRanges object has 3 rows and 6 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.
print(df)
#   Chromosome  Start  End      Name  Score Strand
# 0       chr1    155  200  feature2      0      +
# 1       chr1    900  901  feature4      0      +
# 2       chr1    155  200  feature3      0      -

See the docs for more info: https://biocore-ntnu.github.io/pyranges/


Merge 2 VCFs from different variant callers I am working with WES data for detection of somatic variants and I have used two variant callers because no variant caller is complete in itself. I have used GATK Haplotypecaller for small variants like snp,ins,del etc. And for large variants I have used pindel (specially for SVs detection). Now I have two vcf file for each patients, one from gatk and another from pindel. Although I know various vcf-merging tools but is it right to merge two different vcf file from different callers. Do they give right and proper merged vcf file. Are those vcf-merging tools SV-aware and consider only trivial call for the same variant represented differently in different callers. I know these are too many questions for a single thread but all these are related to each other. I am in dilemma whether I should merge them or use them separately.

Thanks.

The merge utility from SURVIVOR seems to be what you are looking for. I did not try it yet, but it seems to be especially designed to handle SV.

The simplest thing to do would probably to first merge your patients files together for each tool and have one multi-patient file per variant caller. Then I would merge the multi-sample files of the different tools.

Here is a detailed blog post describing the different challenges associated with merging structural variants from different callers and explaining the steps they took to address it. All the code to reproduce the pipeline is available on a github repo.


Subset on multiple genes in Seurat I know that I can do subsetting on just one gene in Seurat:

seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[1], accept.low = 0.1)

However, I want to subset on multiple genes. Is there a way to do that? I just do not want to do manual subsetting on 10 genes, then manually getting @data matrix from each subset, and recreating seurat object afterwards.

I was able to achieve this in the following way:

require(data.table)
cell_names <- vector(mode="character")
for (i in 1:length(neuron_ids)) {
  seurat_subset <- SubsetData(seurat_object, subset.name = neuron_ids[i], accept.low = 0.1)
  genes <- colnames(seurat_subset@data)
  cell_names <- c(cell_names, genes[! genes %chin% cell_names])
}
seurat_subset <- SubsetData(seurat_object, cells.use = cell_names)

Would be interesting to know if Seurat provides such functionality out of the box.


How many models should be generated by rosetta for denovo protein structure prediction? How many models should be generated by rosetta for de novo protein structure prediction of a 382 long sequence? In the tutorial, it says 50,000 models but practically it takes a lot of time (weeks) to generate this number. For instance, I am currently running rosetta on a 382 long sequence and in a week it has generated 1500 models; I am using 16 cpus with Openmpi.

Unfortunately you do in general increase your chance of success by going to the 50k range of models or higher. Projects like Rosetta@home exist because even large research groups can't access as much compute power as they'd like. In this case the long time arises from the combinatorial explosion of possible fragment combinations that need to be explored.

Success is rather protein dependent though - sometimes the right structure falls out after only a few models, other times your fragment library does not cover the protein adequately and you could run Rosetta forever without getting a good structure. Have a look at your 1,500 models and see how similar they look.

382 residues is rather long for Rosetta. Another factor to consider is whether you are using predicted contacts as part of the folding process (e.g. in the PconsFold protocol), which will greatly improve the chances of getting a good model and reduce the number of models you need. In fact the PconsFold group switched to CONFOLD for PconsFold2, I think mainly because it was much faster, so you are not alone in finding Rosetta runs long.


How can I reproduce a manual NCBI search with Biopython Entrez module? I'm trying to make a Biopython script that reproduces a manual search on the NCBI website.

My manual search gives me the following URL:

https://www.ncbi.nlm.nih.gov/nuccore?term=ANOS1[gene%20name]%20AND%20refseq[filter]

And the "Search details" box indicates:ANOS1[gene name] AND refseq[filter].

I get 338 results that I "Send to" "File" as "GI List", save as ANOS1_orthologs_manual.txt.

How can I do the same using Biopython?

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

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

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


cell type specific genes - heatmap using rank based approach I was wondering if some one can help me with pseudo r code for my approach here based on this paper to plot a heatmap
Figure 3

Figure description
The top 40 enriched genes percelltypeare shown in a heat map.Onlyhighlyexpressed genes with FPKM ⬎20 areincludedin this analysis. Fold enrichment iscalculatedasFPKM of one cell type divided by the average FPKM of all other cell types. The majority of these genes showed specific expression by only one cell type, with the exception that some are expressed during more than one maturation stage in the oligodendrocyte lineage.
Similar description in their methods:
Differential expression was calculated as the FPKM of a given cell type
divided by the average FPKM of all other cell types. Genes were ranked by
their fold enrichment in each cell type, and top enriched genes for each
cell type were identified
In my case I have 9 cell types(columns) and ~24000 genes(rows) from mouse and the values are TPMs.
I know to order celltype only based on one column like this
expr.order = expr[order(expr[,1],decreasing =T),]

Which clearly not the right way and would need some help with code in r

The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:

> expr <- matrix(runif(9*24000), nrow=24000, ncol=9)
> top_genes_each_celltype <- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
> dim(top_genes_each_celltype)
[1] 40 9
> top_genes_all_celltypes <- unique(as.vector(top_genes_each_celltype))
> length(top_genes_all_celltypes)
[1] 358

Then you can build a biclustered heatmap with expr[top_genes_all_celltypes,], for example, using heatmap2.


Nextflow error: can not invoke method minus () on null object Hi I was trying to use plink in nextflow and getting the error mentioned in the subject. Here is my code, any help what I am missing. I am new in nextflow, it was ok when input was a single file, but for grouped files, I am getting lost.
params.dir='data/'
params.publish="output/"
Channel
.fromFilePairs("${params.dir}*.{bed,fam,bim}",size:3){file -> file.baseName}
.set {plink-data}

process test1{
publishDir params.publish
input:
set t1,file(plk1) from plink-data
output:
file "${t1}.{bed,fam,bim}" into test1_results
script:
output_file = "${t1}"
base = plk1[0].baseName
"""
module load plink
plink --bfile $t1 --maf 0.05 --make-bed --out ${output_file}
"""
}
process test2{
publishDir params.publish
input:
set file(bed), file(fam), file(bim) from test1_results.collect()
output:
set t2,file '${plk2}.{bed,fam,bim}' into test2_results
script:
output_file = "{t2}"
base = plk2[0].baseName
"""
plink --bfile $t2 --maf 0.01 --make-bed --out ${output_file}
"""
}

Channel names must consist only alphanumeric characters including underscores. Your channel called 'plink-data' contains a hyphen. This is problematic because it is interpreted as a minus() operation to remove 'data' from 'plink'. You get a NullPointerException since 'plink' hasn't been defined. You can test this easily with the following example:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr-names }

chr-names.view()

Results:
Cannot invoke method minus() on null object

The solution is to just use an underscore instead of a hyphen:
Channel
    .of(1..23, 'X', 'Y')
    .set { chr_names }

chr_names.view()

Follow-up issues:
I think if you want to use the bfiles plink option (and I think this is reasonable), I think it makes the most sense to just keep your bfiles in a tuple:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

If you need to reference each bfile separately, you could just unpack the tuple and reference each file accordingly:
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

Note that in the first example, bfiles could be replaced with "*" since no reference is made to the tuple of bfiles anyway. Also note that the list of files produced using the fromFilesPairs operator are returned in lexicographic order.


Subset a multisample VCF file I need some help with GATK. I have a multisample vcf file containing 18125 samples. After performing PCA I have around 303 samples as outliers and I wanted to remove these samples from my multisample vcf file.

I tried doing it using the GATK SelectVariant tool. It does some stuff and generates a VCF file which is a bit smaller in size than the original file, but in the end the total count of samples remains the same. As I have no prior knowledge of GATK (working with it for the first time) I am just wondering if you guys can help me out in fixing this issue.

The command that I used was

gatk SelectVariants -V filtered.vcf -xl-sn remove.txt -O test.vcf

remove.txt contains the sample names in each row (1 sample per row).

After running this command it generated a test.vcf file, for which I counted samples using:

bcftools query -l test.vcf | wc -l

However, I got exactly same number of samples as in the original VCF file filtered.vcf.

Based on the way you've called the command (using -xl-sn), I assume you are using GATK 4 and not GATK 3.

If this is the case, the documentation for the -xl-sn flag is available here: https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_hellbender_tools_walkers_variantutils_SelectVariants.php#--exclude-sample-name.

By default, the -xl-sn SAMPLE means "exclude the sample called SAMPLE from my VCF". In your case, you're providing an entire list of samples, in a file, which is not the default behaviour. In order to tell GATK that you're providing a file and not a sample name, you need to ensure the file you provide has the extension .args.

From the documentation:

  File names must use the extension ".args"


Extracting sequences from FASTA beginning with common 5' end I am trying to figure out the best way to extract sequences from a FASTA file which begin with a common 5' region of 43 nucleotides. Preferably, I would like to to allow for "fuzziness" in this region to allow for mutations or read overlaps. The idea is that any sequences that begin with this region, regardless of size, could be extracted into a new file.

Example of the FASTA file:

>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

And what I would like to get from the code is such:

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

The FASTA file is quite large and I have tried using a couple grep and awk methods retrieve the sequences. Any help you could provide is much appreciated.

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

Save the following python3 code in hit_finder.py

#!/usr/bin/env python3

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

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

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

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

Make it executable: chmod +x hit_finder.py.

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

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

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

>subject
AATCTGTAGTACGACTCACTATAGC

Your query sequences are stored in queries.fa as:

>1-69050-454.08
GTACGGGGAAGGACGTCAATAGTC

>2-65989-433.95
AATCTGTAGTACGACTCACTATAGC

>3-62181-408.91
GATCTGTAATACGACTCACTATAGG

>4-49959-328.53
GGGGAAGGACGTCAATAGTCACAC

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

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

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

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

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

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


Calculating most abundant transcript from RNA-Seq data vcf2maf uses VEP to annotate variants, and I believe selects the default Ensembl transcript to use for annotation. Sometimes the transcript that VEP selects is not the transcript I'm interested in, usually because the selected transcript is not the most highly expressed transcript in my tissue of interest (skin). vcf2maf allows you to provide a transcript override list so that VEP annotates the variant using the specified transcripts instead.
I have several skin samples sequenced with bulk RNA-Seq. I want to estimate the average abundance for each transcript across all samples and then use these abundances to rank transcripts from most to least abundant. Then I will use the most abundant transcript as the default VEP transcript. I plan to use salmon or kallisto to quantify transcript abundance. Should I use TPM or normalized counts to calculate average expression?
My initial thought is to use normalized counts (generated by DESeq2 from raw counts). Are there any problems with this approach? GTEx displays transcript abundance with average TPM, but I thought TPM was inappropriate to use across samples because it doesn't account for between sample differences.
Update: I forgot to mention I also tried using TPM ranks like @ATpoint describes. I haven't fully compared how this compares to transcripts identified by normalized counts, but the initial genes I checked showed good concordance between methods

Transcript abundance quantification is a tricky topic since a read often could belong to several transcripts, so any "count" is a best guess as to which transcript it actually originates from. That being said, there are tools that can help you here:

salmon (as you mentioned) to quanitfy. Run it with --numGibbsSamples 50 (or higher if your computer has the hardware. Also run with -d for the next tool
terminus is a tool that also originates from the COMBINE-lab, It's relatively experimental, but it applies well to your use case. It will take the salmon output and assign counts to transcript groups. If a set of transcripts can be clearly quantified above a certain confidence, that transcript will stay unique. If it is confounded with other transcripts, those will be placed in a transcript group and their counts will be summarised. More info here

You are right in that TPM, FPKM etc. are not appropriate measures to compare across samples. TMM, or DESeq2's median ratio method (what you get with counts(dds, normalized=TRUE) are a step in the right direction, but are still inappropriate across samples. I'd suggest this often posted blog post in the topic.
There are no perfect solutions here since RNA-Seq is inherently a relative quanitifcation method so this will always be an issue. IMO, use DESeq2::varianceStabilizingTransformation(), DESeq2::rlog() or limma::voom() after terminus and rank your transcripts that way.
Just keep the inherent transcript incertainties in mind.


Optimum parameters for genome assembly via canu I have been trying to assemble my genome via canu. Canu gave me about 200 contigs using default parameters. My genomesize is ~9Mbp. Can you recommend some parameters I can tweak to decrease the contigs number?

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

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

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

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

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


probability of an entire sequence, not just one value I am reading "Statistical methods in bioinformatics"

I came across this formula:

As explained in the book O are the observations:

My question is why the author only specified $prob(O)$ for a sequence of length=1 ?, I assume he is only talking about one observation, but generally we have a sequence $O_1 ,O_2,O_3,...$.

I assume that $Prob(sequence)= P(O_1).P(O_2)...$ but I am not sure.

The author writes in the quote that $O$ is to be taken for the whole set of positions $O_{i}$.
In the extracts that you provide there is no reason to believe that it is only for a sequence of length 1. The author writes "the observed sequence $O = O_1, O_2...$"
Thus as you say $Pr(sequence) = Pr(O) = Pr(O_1) * Pr(O_2)...$ according to my understanding of the quotes.


Identifying outlier loci in R I have done some genome scan analysis and I want to identify the outlier loci. For an expected value of the FRD equal to q = 10, a list of candidate loci can be obtained by using the Benjamin-Hochberg procedure as follows:

q=0.1
L = length(adjusted.p.values)
W = which(sort(adjusted.p.values) < q * (1:L) / L)
candidates = order (adjusted.p.values)[w]

and this is part of the outliers and my output

head(candidates)
 [1]  3  4  9   10  13

But I need the output in a different type.
I want to assign TRUE if a locus has been chosen as outlier and FALSE if it has not

  head (desired_output)

  [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE

Any idea how can I modify the Code?

As Devon said, removing the which will give you W as a logical vector that can be used to index the sorted p-values. To get a logical vector that can be used to subset the original (unsorted) values you’ll need to go a step further:

candidates = is.element(1 : L, order(adjusted.p.values)[W])

However, I’m not sure this will work with your data: the name adjusted.p.values suggests that these values are already FDR-adjusted. So you cannot/should not perform additional correction on them. If they are already adjusted, then you can get candidate loci simply by

candidates = adjusted.p.values < q


How to extract RNA sequence and secondary structure restrains from a PDB file I'm trying to find a programmatic way to automatically extract the following information from a PDB file:

RNA sequence
Secondary structure restraints in bracket format, e.g. . (( . ( . ) . ))

Does software exist that can take a PDB file as input and generate these two pieces of information?

e.g. file. 3NDB_ba.pdb

The nucleotide sequence is:

gUCUCGUCCCGUGGGGCUCGGCGGUGGGGGAGCAUCUCCUGUAGGGGAGAUGUAACCCCCUUUACCUGCCGAACCCCGCCAGGCCCGGAAGGGAGCAACGGUAGGCAGGACGUCGGCGCUCACGGGGGUGCGGGAC.

And the secondary structure :

.(((.(..(((((((((.(((((..(((((.(((((((((....)))))))))..)))))....((((((.....(((.....(((....))).....)))..)))))).))))))).))))))).....).))).

I suggest you take a look at rna-pdb-tools we do way more than you need! :-) The tools can get you a sequence, secondary structure and much more using various algorithms, and all is well documented http://rna-pdb-tools.readthedocs.io/en/latest/

To get sequence http://rna-pdb-tools.readthedocs.io/en/latest/main.html#get-sequence

$ rna_pdb_tools.py --get_seq 5_solution_1.pdb
> 5_solution_1.pdb A:1-576
CAUCCGGUAUCCCAAGACAAUCUCGGGUUGGGUUGGGAAGUAUCAUGGCUAAUCACCAUGAUGCAAUCGGGUUGAACACUUAAUUGGGUUAAAACGGUGGGGGACGAUCCCGUAACAUCCGUCCUAACGGCGACAGACUGCACGGCCCUGCCUCAGGUGUGUCCAAUGAACAGUCGUUCCGAAAGGAAG

or you can get sequence and secondary structure via x3dna (http://x3dna.org/)

[mm] py3dna$ git:(master) ✗ ./rna_x3dna.py test_data/1xjr.pdb
test_data/1xjr.pdb
 >1xjr nts=47 [1xjr] -- secondary structure derived by DSSR
 gGAGUUCACCGAGGCCACGCGGAGUACGAUCGAGGGUACAGUGAAUU
 ..(((((((...((((.((((.....))..))..))).).)))))))

The problem is not trivial. I hope rna-pdb-tools has stuff that can solve this issue for you (http://rna-pdb-tools.readthedocs.io/en/latest/want.html#module-rna_pdb_tools.utils.rna_x3dna.rna_x3dna)


Error while using E-utilities on slurm HPC I keep getting the error while using E-utilities on Slurm HPC, could you suggest how I can resolve it?
SCRIPT:
#!/bin/bash
#SBATCH --job-name=tax_job_test      # Job name
#SBATCH --mail-type=END,FAIL         # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=na@gmail.com    # Where to send mail
#SBATCH --ntasks=8                  # Number of MPI ranks
#SBATCH --nodes=1                    # Number of nodes
#SBATCH --time=23:05:00              # Time limit hrs:min:sec
#SBATCH --partition=k2-medpri

module load e-utilities/03.02.19
module load apps/perl/5.20.2/gcc-4.8.5

for i in `cat taxonomy.id.txt` ; do efetch -db nuccore -format gbc -id ${i} |xtract  -pattern INSDSeq -element INSDSeq_accession-version INSDSeq_taxonomy INSDSeq_organism >>taxo.2020.txt; done

input taxonomy.id.txt is like this:
AKN43077.1
AKN41814.1
PNE10086.1
AKN42672.1
AKN41654.2
AKN43077.1
PNE10086.1
AKN43189.1
TMK08501.1
TMK08501.1
AKN42976.1
AKN43189.1
AIZ47938.1

Error:
The script keep giving below error:
501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)
No do_post output returned from 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network'

Result of do_post http request is

$VAR1 = bless( {
                 '_rc' => 501,
                 '_request' => bless( {
                                        '_method' => 'POST',
                                        '_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
                                        '_content' => 'db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network',
                                        '_headers' => bless( {
                                                               'content-type' => 'application/x-www-form-urlencoded',
                                                               'user-agent' => 'libwww-perl/6.39'
                                                             }, 'HTTP::Headers' )
                                      }, 'HTTP::Request' ),
                 '_content' => 'LWP will support https URLs if the LWP::Protocol::https module
is installed.
',
                 '_headers' => bless( {
                                        'client-warning' => 'Internal response',
                                        '::std_case' => {
                                                          'client-warning' => 'Client-Warning',
                                                          'client-date' => 'Client-Date'
                                                        },
                                        'client-date' => 'Mon, 05 Oct 2020 14:27:37 GMT',
                                        'content-type' => 'text/plain'
                                      }, 'HTTP::Headers' ),

Many thanks

501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed) indicates that a specific perl module is missing.
If indeed your edirect installation is from March 2019 (as I interpret from e-utilities/03.02.19), you should update your edirect installation to the latest version. If you update, I expect the necessary perl module will be installed and the error will go away.


Read length distribution from FASTA file I have a ~10GB FASTA file generated from an Oxford Nanopore Technologies' MinION run. How can I quickly and efficiently calculate the distribution of read lengths?

A naive approach would be to read the FASTA file in Biopython, check the length of each sequence, store the lengths in a numpy array and plot the results using matplotlib, but this seems like reinventing the wheel.

If you want something quick and dirty you could rapidly index the FASTA with samtools faidx and then put the lengths column through R (other languages are available) on the command line.

samtools faidx $fasta
cut -f2 $fasta.fai | Rscript -e 'data <- as.numeric (readLines ("stdin")); summary(data); hist(data)'

This outputs a statistical summary, and creates a PDF in the current directory called Rplots.pdf, containing a histogram.


If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it? I am reading a book about RNAseq analysis and it says

"To calculate the probability that a read will map to a specific gene, we can assume an average gene size of 4000 nt (100 M nt divided by 25,000 genes). At 30 M reads equivalent to 30× coverage, at single read 100 nt (or paired-end read 50 nt) length, we can expect a single read to map to the average expressed and length gene 4000 nt× 30 coverage/100 nt 1200 times. Thus, if the gene is expressed at a level of 1/1200 compared to the average gene, then we have a 50:50 probability to have a read map to it."

I understand that a single read would map to the average expressed 1200 times but I don't understand how the 50:50 probability is calculated. Please could you explain how this probability is determined.

I actually disagree with that.. I guess I write this down as a discussion.
The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read.
However, because of sampling, sometimes you get 0, 1 or 2. If we assume the sampling follows a poisson distribution, we can calculate the probability of getting 0, 1 , 2 ... reads:
plot(dpois(0:10,lambda=1),ylab="Probability",xlab="No of reads")

The probabilities are like:
data.frame(n=0:5,p=dpois(0:5,lambda=1))

  n           p
1 0 0.367879441
2 1 0.367879441
3 2 0.183939721
4 3 0.061313240
5 4 0.015328310
6 5 0.003065662

The probability of getting exactly 0 or 1 read should be ~0.37. The probability of getting at least 1 read would be 1 - P(n=0) = 1- 0.367879441 = 0.6321206


Pandas automatically rounds GWAS P-value I am working with a specific GWAS.  If I were to run this on the command line

grep <rs_id_of_interest> GWAS.txt

I would see the GWAS p-value to be on the order of 10^-200. Very tiny p-value.

However, if I were to do the same but in pandas, with the following code,

gwas = pd.read_csv(GWAS_file,sep='\t')
gwas_rs_id = gwas.loc[gwas['marker'] == 'rs_id'].copy()

and then print gwas_rs_id, it shows the gwas p-value to be 0.0. Is there a way to read the file and not have it automatically round the tiny tiny p-values?

Thanks!

The solution you likely want is here,

pd.set_eng_float_format(accuracy=x, use_eng_prefix=True)

x = whatever is required
The function set_eng_float_format has been moved around a bit and is now a top level function

You might be dealing with maximum likelihood, so convert to log likelihood -200 that sort of thing

Best idea

gwas_rs_id['logged'] = np.log(gwas_rs_id.marker)

First idea,

import decimal as D
....
gwas_rs_id['marker'] = gwas_rs_id['marker'].astype(D.decimal)

Last attempt,

pd.set_option('display.precision', x) # x whatever dp you need, this might work in combination with a log likelihood

Try a gwas_rs_id.dtype() check whether you're using float32, float64 or still in object. It should convert to scientific notion (E values) automatically.


RAD Seq Data Analysis without barcode I received a ddRAD Seq data from my supervisor without barcodes and restriction enzymes. I asked him for both and he said I don't need it since the data has been cleaned by the company. Now, I am kinda confused and don't know where to start.

Wow Biobash, that's kinda harsh. It is cool that sometimes you unexpectedly learned something new. What Michael mentioned is pretty interested, no one has ever mentioned that or I have never heard that.
So, my approach would be that you now are doing a Poolseq study, since you can't identify individuals. The other way to try to catch the paralogues is to map your reads to a sister species. My question is what are you trying to estimate? Genetic differentation? Genetic diversity? Pop Gen structure?
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6116966/
https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.5240
https://www.g3journal.org/content/9/12/4159
One semi-robust approach to assess usefulness is thinking about your system and defining what you expect to find, let's say in terms of genetic differentiation or genetic structure, and estimating that from your data. If that hypothesis is confirmed, your dataset might be useful. You can test the same by first estimating the loci by mapping your reads to a sister reference genome, to try to detect paralogues. Since you don't have barcodes, then each "locus" is independent, and you lose power, but there is nothing else you can do.
Good luck;


Enrich for genes expressed in the plasma membrane I have a list of cell markers differentially expressed between 2 groups of cells, that I got with the FindMarkers command in Seurat. I want to enrich for those genes that are expressed in the cell membrane. Are there any specific tools/commands I can use to do that? Thanks!

You can do functional enrichment analysis, for instance with gProfiler in R and see if your marker genes are enriched in any GO term related with the plasma membrane.

Alternatively, you can generate a list with genes whose product are located at the plasma membrane and perform Gene Set Enrichment Analysis, for instance with fgsea in R.


Counting reads within the intron I have paired-end RNA-seq samples and a list of intron coordinates in bed format. I want to count the intronic reads such that:

1. The reads should overlap the intron by atleast 25 bp
2. If both the reads of the pair are within the intron, count it as 1.

What tool can I use to compute it.

Thanks

You could use BEDOPS bedmap to map reads to introns, using 1) the --count operator to do counting of reads overlapping by your criteria; and, 2) the --indicator operator to do a true/false operation, where reads are contained entirely within the intron.

For instance, to count reads that overlap introns by at least 25 bases, use --count to count and --bp-over 25 to specify the minimum overlap threshold:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

To filter for introns where reads are contained entirely within the intron, use --indicator to get introns that meet the criteria of full containment of reads, via --fraction-map 1:

$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- > answer.indicator.bed

Once you have answer.counts.bed and answer.indicator.bed, you can compare those two results with bedmap --indicator --exact to see where there are commonalities (true/1) and differences (false/0), and increment the final count where there are identical introns, e.g.:

$ bedmap --indicator --echo --exact --delim '\t' answer.counts.bed answer.indicator.bed | awk '{ $NF += $1; print $0; }' | cut -f2- > answer.final.bed

If your reads are in a different format, i.e. BAM, you can use bam2bed in a process substitution, replacing reads.bed with <(bam2bed < reads.bam), e.g.:

$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed <(bam2bed < reads.bam) > answer.counts.bed

Etc.

You can also just convert separately, so that you only need to do it once:

$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

Etc.


How to extract fasta from a blastdb How to extract the sequence used to create a blast database.  This is useful when you download a blastdb from somewhere else e.g. one of the databases provided by NCBI including the 16SMicrobial database. Or alternatively, when you want to double check which version of a sequence you have included in a blastdb.

You can extract fasta sequence from a blastdb constructed from a fasta file using blastdbcmd which should be installed when you install blast/makeblastdb.

blastdbcmd -entry all -db <database label> -out <outfile>

If you had a database called my_database which contained the files:

my_database.nhr
my_database.nsq
my_database.nin

and you wanted your fasta output file to be called reference.fasta you would run the following:

blastdbcmd -entry all -db my_database -out reference.fasta


How do I get repeating regions from fasta file or directly from internet If I go to ncbi site, apply customize view -> Customize -> All features -> Update view and then ctrl + f : Alu I can get all alu repeaing regions coordinates from this chromosome.
It looks like:
 repeat_region   complement(37490..37800)
                 /note="AluY"
                 /rpt_family="SINE/Alu"
                 /rpt_type=DISPERSED

How could I download all this coordinates and names of repeating regions from the internet to .txt file? Or some tool

Click on the 'GenBank' dropdown menu at the top left corner of this page to see a graphical view of the sequence. You should see a 'Repeat region' track in the Sequence Viewer. Click on the 'Download' button in the top right corner of the Sequence Viewer, choose 'Download Track Data' to launch a dialog box where you can enter a range, select the 'Repeat region' track and download data in BED format. Unfortunately, the names for the features are not great.
If you want the names to be the same ones in the flat file, you can download the feature table and parse it. To download the feature table, click on the 'Send To' link at the top right corner of this page, choose 'Complete Record', 'File' and 'Feature Table' as 'Format'. After you download the feature table, you can parse it using awk as shown below to create a 4-column table with the data.
cat seq.ft \
  | awk 'BEGIN{FS="\t";OFS="\t"}{
    if(NF==3) {print s,e,f,q ; s=$1; e=$2; f=$3; q=""}
    else if(NF==2) {print s,e,f,q ; s=$1; e=$2; q=""}
    else if(NF==5) {q=q"|"$4"="$5} }' \
  > feat_table.tsv


RNA seq .counts.txt to bigwig conversion I've downloaded neuronal RNA seq data from GEO. The files are in .counts.txt format. Could I convert them to bigwig? If so, how?

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


Model the effect of treatment on mutation burden I am struggling a bit to model the following problem. Basically, I would like to model tumor mutation burden (mutations per megabase, a continuous) as a function of treatment (categorical). The following is an example of how data looks like

patient-id therapy-1 therapy-2 therapy-3 cancer-type TMB
PID-1      1         NA        1         CT-1        1.3
PID-2      0         0         1         CT-1        0.14
PID-3      NA        0         1         CT-2        5
PID-4      1.        1         NA        CT-3        2.3
PID-5      1         1         NA        CT-3        10
...
...

1 - treated; 0 - not treated; NA - No information

I was trying multiple linear regression, but in most cases one patient got multiple therapies. This is where I am confused on how to interpret the results of linear regression. So far, I haven't included the cancer-type information into the model. However, IMO, including cancer-type into model will yield specific treatment effects in specific cancer type? Or is it over-expectation of what we can achieve from linear model?

Also, is there a better way to address the problem to get reasonable associations between treatment and TMB in this case?

Thank you for the help. Stay safe!

Depends how you consider "linear regression", if you mean Pearsons, Spearman's analysis then general linear modelling (GLM) is a far better approach if you have continuous data. I don't really understand TMB, but that is a critical output then GLM will work. Logistic regression is an alternative approach data type depending. You have to perform a Q-Q plot and the residual against against any transformations you make (you might need to transform TMB). The lower the residual the better, the straight the Q-Q plot the better.

To explain how GLM works long hair and height correlate in a straight regression analysis. The implication is if you want to be tall - cut your hair! Of once you account for gender and send it through a GLM then gender and hairlength will correlate far stronger than hair length and height.

Mac OSX has a very easy to use pointy clicky, but highly professional, GLM package on its AppStore called Wizzard (bit expensive). It doesn't do general linear mixed modelling however. Everything does GLM of course, its just this package is really easy to use and you can figure out GLM by trial and error.


Tree cut height for Consensus data input not working I have tried different combination of tree cut heights but im getting the exprSize as one sample as full the other one as empty .

# Choose the "base" cut height for the female data set
baseHeight = 135
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(135, 145*exprSize$nSamples[2]/exprSize$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
  plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
       xlab="", sub="", cex = 0.7);
  abline(h=cutHeights[set], col = "red");
}
dev.off();

My figure

When i see the cutHeights and exprSize i do see its higher than 120 but still i get only one set of sample which is 55 whereas the my other set which is

exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55 47

$structureOK
[1] TRUE

This is after tree cut

exprSize
$nSets
[1] 2

$nGenes
[1] 8213

$nSamples
[1] 55  0

$structureOK
[1] TRUE

> cutHeights
[1] 135.0000 123.9091

Im not able to figure out what is wrong

Your data don't have obvious outliers so you don't have to worry about removing them. The problem is that both sets show prominent sample clusters. I would try to figure out what drives those clusters (is it technical or biological), and whether the variation should be removed. See the comments in WGCNA FAQ under point 5 (My data are heterogeneous. Can I still use WGCNA?). This looks like data from GEO, so look into the sample annotation for these data sets and see if any of the recorded sample characteristics align with the clusters you see in the sample trees.


Filter reads belonging to unique sequences with threshold I have CLL samples as fastq files and I want to remove those reads which have a unique sequence with less than 10 read counts each. I tried this by following way using awk to make it faster:

for fq in $(ls fastq1/*.fq.gz); do
    basename=$(basename "$fq" | cut -f1 -d '.')
    zcat $fq | awk "/^@/{getline; print}" |sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt; grep -B 1 -A 2 -f tmp.txt $fq | grep -v ^--'| pigz -p 10 > fastq2/$basename.fastq.gz

done

rm tmp.txt

But this is not working well as many reads with more than 10 read counts are also missing. I think it's OK to stick with bash/awk for fast fastq file processing. What am I missing here, or is this command wrong. Can anyone suggest me the right and fast way to do it??

Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.

I'd adjust your script to the following:

for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
 zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

Instead of

for fq in $(ls fastq1/*.fq.gz)

use this

for fq in fastq1/*.fastq.gz

Iterating over ls is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.
Rather than removing the separators that grep will introduce by filtering them out, use the following option to not create them:
grep --no-group-separator
And, minor bug, as others have mentioned, you should use zgrep or zcat | grep since your files are compressed.

Oh - and quote your variables.

Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.