Where can I visualize Ramachandran plot? I am using Rampage webserver to obtain a Ramachandran plot. However, the server does not reveal any plot just values (favored, allowed, outlier regions).

Can you recommend other similar web sites?

Here is a page with several suggested methods. I agree that Rampage looks the best of the options, but it seems that it provides some suggested workflows for drawing it in R so that it is very customizable and you can probably get it to show any features you want if you're willing to customize a little.

Here is a Python tool for drawing Ramachandran plots.

There are additionally several other servers/tools.


How to resolve in snakemake error : "Target rules may not contain wildcards." I would like to do easily reproducible analysis using publicly available data from NCBI, so I have chosen a snakemake.

I would like to write a single rule, that would be able to download any genome given a species code name and separated table of species and their NCBI IDs. So I wrote a script scripts/download.sh that takes a species code and download the genome to data/<species_code>/genome.fa.gz. The script internally reads the table tables/download_table.tsv, where are corresponding species code names and NCBI IDs.

So I tried to do a snakemake like this :

species='Cbir Avag Fcan Lcla Dcor Dpac Pdav Psp62 Psp79 Minc1 Minc2 Mjav Mare Mflo Mhap Pant'

rule download:
    input :
        "tables/download_table.tsv"
    output :
        "data/{sp}/genome.fa.gz"
    shell :
        "scripts/download.sh {sp}"

However, snakemake returned an error message I do not really understand :

Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards.`.

Is there a way to write a single rule for downloading all the genomes?

The problem is that you need a master rule that requires all of your desired outputs as inputs, in your case it would be :

rule all:
    input:
        expand("data/{sp}/genome.fa.gz", sp=species.split(' '))

You'll also need separate download link inputs for each species. You could make a separate download_table.tsv for each species, but it would probably be easier to make a config file with this information, and add a params keyword to your rule. Something like:

rule download:
   params:
       url=config['locations']['sp']
   output :
       "data/{sp}/genome.fa.gz"
   shell :
       "scripts/download.sh {params.url}"


Metapath2vec on Drug-ADR Heterogeneous Graph I am new in this field and I am having some problems regarding a new project.
I built a graph using Drugbank Data connected to SIDER Adverse Reactions.
I used Organ- level Terms to classify the ADRs in 25 different groups.
I aim to use the algorithm Metapath2vec in order to cluster drugs based on the group of ADRs they lead to.
I am currently using the python package for Machine Learning on graphs Stellargraph and the related implementation of Metapath2vec.
I chose Clustering as downstream task with the resulting node embeddings, e.g. DBSCAN but the results are not promising.
I would like to have some clusters of drugs related to different adverse drug reaction.
Since I am new in this field everything I am trying is based on scientific literature but I don't know if this is the right approach for my objective.
This is the code related to the metapath2vec algorithm:
walk_length = 100  # maximum length of a random walk to use throughout this notebook

# specify the metapath schemas as a list of lists of node types.
metapaths = [
    ["drug", "adr", "drug"],
    ["drug", "adr", "drug", "drug"],
    ["drug", "drug"],
    ["drug", "adr", "group_adr", "adr", "drug"],
    ]

# Create the random walker
rw = UniformRandomMetaPathWalk(graph)

walks = rw.run(
    nodes=list(graph.nodes()),  # root nodes
    length=walk_length,  # maximum length of a random walk
    n=1,  # number of random walks per root node
    metapaths=metapaths,  # the metapaths
)

from gensim.models import Word2Vec

model = Word2Vec(walks, size=128, window=5, min_count=0, sg=1, workers=2, iter=1)

Should I change approach?
what other algorithm for representation learning can I use in order to reach my goal?
what can I improve in the presented approach to have better node embeddings?

Clearly the additional references do put Metapath2vec into good use. Its a many to many ML relationship. Most catagorization/classification particularly in deep learning is many  to one. So thats the good thing about the technique you  are using, you've 25 catagories, which is quite large for machine learning(ML) and large for deep learning. The technique appears to be half way between supervised and unsupervised learning. The random walk establishes a zone (circle) of catagorisation.
NLP Background
The situation where it would work in NLP (natural language processing) is for example a spell checker, where you have one correct word and loads of possible typos, some of which you not observed in training. So what will happen is the random walk will pick up the local vacinity around the 'correct spelling' by the vectorization (converting a word into numbers) of the incorrect spelling will land it close to 'true spelling' in a 2D vectorization plot. The training following the random walk says  'yep you belong with me'.
The papers given talk about being able to identify the context of a word within a sentence which is a very NLP thing, one vectorisation being 'bags of words'. So here the sentence represents one of many categories of information, i.e. each is the central node within the graph, the vectorisation of the sentence places of meaning closer to one 'hub' than another and again the local vacinity is identified by the random walk. In training this 'random walk vacinity' is identified, so it will algorithmically recreate the 'random walk locality' so that sentences that have not been trained can have their meaning identified.
The central theory of deep learning, which not quite true of other ML, is any algorithm can be recreated through training. So like it will sort of recreate a 'random walk zone'
Difference with drug discovery
The papers you have map drugs to genes, but it means you have to have a "true drug to gene" relationship. The problem I see here is what is 'true' for a drug-protein interaction? With a spell-checker there is one 'true spelling', with a sentence there is pretty much one 'type of meaning', like a movie was great, average, terrible based on a ratings score and stuff like gender classification. If its not 'true' the random walk will get it wrong. So its a big assumption. If you have a circle and pick a point on the edge of the circle then try and draw a circle around the original circle, you can't unless the circle is so large it risks including stuff that is not part of the original circle. The problem the papers is how do they know what is the 'centre' of the circle? Its an a priori assumption and you need really good biological rationale.
With your stuff I agree its easier because adverse drug reaction, some drugs could be fatal or toxic (traditional cancer drugs or parasite drugs), so you'd pick the worst case adverse reaction. That then centres the definition around which everything else is relative. If however your 'true catagory' turned out to be an outlier its not going to work, because the random walk only works if its starts bang in the centre of the 'catagory'. Random walks form spheres, in case you are not aware. I can see this working as an advanced 'unsupervised learning' technique.
One very real advantage is if a drug doesn't have an adverse reaction it MIGHT (I'm guessing) not give a result, in which case you are good to go - that a good drug. If that is what happened that would be amazing because its not how ML, or particularly DL works. Traditional ML/DL will force a classification even if it doesn't fit. There are methods to accept some sort of "half-classification" but I don't know how to do it. This is one of the reaons you need lots of things to train.
I am assuming that vectorisation of your drug is fine because there is a huge history of this in drug-discovery going right back to the origins of the field.
Summary
So yeah I get the idea, I see the flaws and the advantages. The final issue is that most ML requires pretty large data sets and DL require huge data sets and when you've got 25 catagories. For DL approaches the data set size might need to be huge unless you can do 'transfer learning' (backing your training onto a pre-trained model). If your method resulted in a zero classification i.e. "doesn't fit anything I've got" - that would be cool and that is the key advantage I can see in the 'random walk' approach which is not available using traditional methods.
As long as you are aware, you do not have to use this algorithm to get a classification. This final point what appears missing in the papers is using a 'traditional' method by comparison. The ML before deep learning, stuff like random forests, routinely takes loads of catagories. If you are using a non-deep learning ML you would traditional use several different approaches to justify your result.

My answer is a bit rambly but I think you'll get the gist. Peronally I would complement with a standard ML, but all methods have their flaws. Deep learning (DL) has a big history in drug discovery (alot of work was done via Tensorflow), but I'm not sure you've data set size for this.


Finding drugs based on gene expression data I've got a list of genes and I want to assess their association with drugs based on gene expression data sets.

A comprehensive tutorial would also be helpful.

Any help would be appreciated.

It would be great if you can provide more information about the gene expression data set. If you have response to single drug, the general approach is to find the deferentially expressed genes with tools like DESeq2, or make a heat map based on the expression of your listed genes to see the drug effect before and after the treatment. But, if you want to study the effects of combination of drugs it is necessary to consider the additive and synergistic effect of drug combinations.


Tracing and quantifying inheritance over many generations Apologies if this is the wrong location for this question, and please feel free to point me in the right direction with appropriate Anglo-Saxon embellishments...
I'm after some ball-park numbers for the following questions (plus any hints how to calculate them for different params - e.g. less generations, larger starting pop.):
Given a starting population of 10,000, and given 2,500 generations, with an end population of 5,000,000,000

What are the chances of being a direct descendant of 1 individual of the original 10,000?
If you had a DNA sample of that individual, would you be able to tell if you were a direct descendant?
Would siblings of the original individual confuse things with regard to Q2?

Since it seems my comment helped, I've been suggested to post this as an answer...
I'd reccomend everyone interested in these sort of population genetics to have a look at 'The Ancestors Tale' by Richard Dawkins. In particular, chapter 0, which deals with human/hominid genetics. That should help answer questions. The relevant part of the book appears to be available online for free in GoogleBooks (https://books.google.co.uk/books/about/The_Ancestor_s_Tale.html?id=s86_CQAAQBAJ).
Below I'll note down the answers from the comments section, not because I wrote them but I figured it might help anyone who stumbled over this question afterwards:

The chances of being a direct descendant of 1 individual of the original 10,000 is 100%, assuming they have any descendants. - Tom Melly

If one had a DNA sample of that individual, one could use uniparental markers (either Y chromosome for men, or mtDNA for women) to infer whether they are a direct descendant. - 4galaxy7

(My guess is that it doesn't have to be uniparental markers. Other genes might work just as well. I'm thinking of Neanderthal/Denisovan work. But I might be wrong.)


Find a cutoff value for genes that are expressed in single cell RNA-seq? I want to find a cutoff value for each gene, above which we can consider a gene expressed.

The problem is that not all effectively non-expressed genes will have 0 counts due to sequencing errors for example.

My question is related to this one, but is concerned with single-cell data. The technical differences between bulk and single cell RNA-seq warrant a separate question for this problem. Also similar to this one, but focuses on gene expression and not differences between zero counts.

Not sure if this has anything to do with your question but I think the way how Monocle approach this question sort of makes sense to me. The normalization methods for Seurat and Monocle are somehow similar (correct if I am wrong, this is simply my personal experience), as @kohlkopf mentioned.

For monocle specifically, it has this function where you can display dispersion v.s. expression of the genes to decide the cut off. Here is an example of this plot that I generated when processing the PBMC data downloaded from the Seurat tutorial.

You can assign the cut-off manually but I doubt I answered your question on how to decide the cut-off value. Maybe you need to try a set the expression cut-off from a range to experiment what value may potentially answer a biologically relevant question in your context.

If you are interested in how monocle approach this, I think the Monocle introduction and Dave Tang's Blog are excellent in explain how it works.

Again, you are right that 0 could be sequencing error. But I would also suggest that, if you are looking at genes at very small expression level, scRNA-seq, at least droplet-based technique, could be really difficult in such context. I happened to hear someone trying to find lncRNA expression talking about how a "0" sometimes may be a true "0" due the low expression nature of lncRNA in scRNA-seq, etc.

Anyhow, not sure if I answer the question but I do hope it helps you somehow. I am not able to make comments but this is the best I can do.


How to introduce point mutations in a protein structure? I want to introduce point mutations to a protein PDB file and automatically predict the 3D coordinates of a new amino acid residue in the mutated protein.
I used interactive Schrodinger's Maestro tool to mutate residues and minimize them to avoid clashes, however this time i want to automate the task and use some open source tools.
Is there any free software for this task?

PDBFixer and pyMutate are simple, lightweight, free and open-source. Both don't do any optimizing, though, so you may have to run some molecular dynamics afterward to get what would be the right "orientation".
VMD would work for sure, I believe PyMOL too. They are both graphical interfaces but you can start them in text only mode (or just import pymol in a python script and use the Python API). They are free for particulars and academics, but given their renown you may already be using them.
If you're looking for something a tad more complete and that may adapt the side-chain to the surrounding environment I'd go with Modeller, which provides two example scripts to do point mutations (here and there).
If you don't want to just get the side-chain of this amino-acid right but actually strive to see some general backbone adaptation too then use PyRosetta. It can probably do what all the others do, too, but given it is 1 GB big it may be overkill. Any of the package coupled with molecular dynamics would also do the trick.


Library for computing BWT-based alignments I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. As far as I can tell, though, both of these tools and their derivatives are invoked strictly via a command-line interface. Are there any libraries that implement BWT-based read alignment with a C/C++ API?

Python bindings would also be great, but probably too much to expect.

First, let us remark that there exist several hundred read mappers, most of which have been even published (see, e.g., pages 25-29 of this thesis). Developing a new mapper probably makes sense only as a programming exercise. Whereas developing a quick proof-of-concept read mapper is usually easy, turning it into a real competitor of existing and well-tuned mappers can last for years.

It is not clear from the provided description how long is your reference, how many alignments you need to compute, etc. In certain situations, it may be useful to write a wrapper over existing mappers (e.g., using subprocess.Popen for running the mapper and PySam for parsing the output); while in some other situations, standard dynamic programming may be sufficient (e.g., using SSW with its Python binding).

Let assume that you want to develop a toy read mapper. Most of read mappers are based on a so called seed-and-extend paradigm. Simply speaking, first you detect candidates for alignments, usually as exact matches between a read and the reference (using either a hash table or some full-text index – e.g., BWT-index). Then you would need to compute alignments for these candidates, typically using some algorithm based on dynamic programming, and report the best obtained alignments (e.g., the ones with the highest alignment score).

There exist two big, powerful and well debugged libraries implementing BWT-indexes which can be easily used for building a read mapper:

SeqAn. See the FMIndex tutorial and the Pairwise Sequence Alignment tutorial for quick examples of how to detect exact matches and how to do pairwise alignments. Also, they provide a tutorial about a quick development of a read mapper, but the resulting mapper seems to be hash-table based, not BWT-based. SeqAn is used, e.g., in YARA.
SDSL-Lite. This is a general library for succinct data structures. See the tutorial slides for an idea of how it works. For instance, GramTools are based on SDSL.


Where to get '--fldMean' and '--fldSD' for single-end Salmon run I need to run single-end bulk-RNAseq salmon - based alignment and I am not sure what I should supply as --fldMean and --fldSD which seem to be crucial for the single-end experiment alignment. Can they be somehow derived from the FASTQ files that I have, or should I be provided with them by someone who did the experiment?

From looking at the salmon documentation here, it is correct that both of these are essential for single end experiments:

  Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads

From the discussion here it would appear that you can't estimate these parameters at all via the FASTQ files alone.

The fragment length refers to the physical molecule of D/RNA captured in the library prep stage which you have sequenced (in part) on your sequencer.  Thus ideally this sort of data should come from a bioanalyzer trace from whoever performed the library prep or similar QC step, this should ideally show the mean fragment length and standard deviation.


How to deal with heterozygosity during polishing of genome assembly based on long reads? All the long-read sequencing platforms are based on single-molecule sequencing which causes higher per-base error rates. For this reason a polishing step was added to genome assembly pipelines - mapping raw reads back to assembly and correcting details of the assembly.

I have decent PacBio RSII dataset of single individual genome of heavily heterozygous non-model species. Assembly went well, but when I tried to polish the assembly using quiver it could not converge over a couple of iterations and I bet it is because of too great divergence of haplotypes.

Is there any other way to polish a genome with such properties?
For instance, is there a way to separate long reads by haplotype, so I could polish using one haplotype only?

A few possibilities:

Falcon

Try falcon and falcon-unzip. These are designed exactly for your problem and your data: https://github.com/PacificBiosciences/FALCON

Not Falcon

If you think you have assembled haplotypes (which seems reasonable to expect given enough coverage), you should be able to see the two haplotypes by just doing all pairwise alignments of your contigs. Haplotypes should show up as pairs of contigs that are MUCH more similar (even with a lot of between-haplotype divergence) than other pairs. Once you have all such pairs, you can simply select one of each pair to polish.


how to easily make complex Euler/Venn diagrams in R? From my microarray data I have selected three subsets of probes. I would like to display these subsets as a Euler diagram to visualize the overlaps, but am struggling with inputting the combinations into Eulerr. All three subsets need to be nested within the full set of probes; one subset is fully contained within another, and the third overlaps with the other two. Is there a package into which I could feed the vectors of probe names and let it determine the overlaps for me? Many thanks in advance!

Since you are in R, consider using the UpsetR library.  It doesn't make Venn diagrams but it helps to visualize overlaps between any number of groups.
http://gehlenborglab.org/research/projects/upsetr/#:~:text=UpSetR%20is%20an%20R%20package,based%20on%20groupings%20and%20queries.


How do I filter a GFF file by gene type? I have a file with following information:

chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    141474  173862  .       -       .       gene_id "ENSG00000241860.6"; gene_type "processed_transcript"; gene_name "AL627309.5"; level 2; tag "ncRNA_host"; havana_gene "OTTHUMG00000002480.3";
chr1    ENSEMBL gene    157784  157887  .       -       .       gene_id "ENSG00000222623.1"; gene_type "snRNA"; gene_name "RNU6-1100P"; level 3;
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

In the above file there is a gene_type with different names like processed_pseudogene, lincRNA, snRNA, unprocessed_pseudogene, antisense_RNA, processed_transcript.

I want to create a new file based on gene_type. Among all the gene_types I want the new file to be with certain names of gene_type like only processed_pseudogene, lincRNA, unprocessed_pseudogene, antisense_RNA.

The new file should look like below:

chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

grep -E 'gene_type "(processed_pseudogene|lincRNA|unprocessed_pseudogene|antisense_RNA)";' input_file.txt

Explanation

grep -E: use extended regular expression so you don't need to escape things.
This kind of expression (foo|bar|baz) allows grep to include either of these strings. Hence, the total expression will match gene_type "processed_pseudogene";', gene_type "lincRNA";', etc.


How to generate Beta diversity boxplot from phyloseq object? I am working with profiled metagenomic taxonomy abundance data. I want to generate beta diversity (bray-curtis) boxplot from a phyloseq object where two groups (control and test) will be shown. Something like this figure:

How can I make it? Can anyone share any code or tutorial?
Thanks,
dpc

You could convert the phyloseq object to a dataframe and just plot it using base R or ggplot2. I'm not exactly sure what the bray curtis df will be called but you can access it using the @ notation then convert with data.frame e.g.
bray_for_plot <- data.frame(phylo_obj@bray)


What is the difference between a MeSH concept and a MeSH descriptor Medical Subject Headings (MeSH) is a very comprehensive vocabulary to index the biomedical literature on Medline or PubMed (which is the major Biomedical literature database). MeSH is arranged into Descriptor > Concept > term hierarchy, but I still do not understand the difference between a MeSH descriptor and a MeSH concept.

For example, for Knee

MeSH descriptor is D007717
MeSH concept is M0012066 (Under the concepts section in the link)

Even after reading the description here and here, I am still not sure about the purpose of MeSH concept and MeSH descriptor.

Could anyone please shed some light on it?

MeSH is organized into Descriptor > Concept > Term hierarchy, where all the synonymous <term>'s are grouped together under a <concept>. For example, consider the MeSH terms,

Blood-Brain Barrier
Hemato-Encephalic Barrier
Brain-Blood Barrier

All these terms are group under the umbrella concept of "Blood-brain barrier".


Function to find ATP binding sites I'm currently doing a python project where I need to write code that extracts the ATP binding site from a FASTA file using a regular expression.

I've tried the following code but it does not work. It opens a new file but the file is empty. I don't know what's wrong.

from Bio import SeqIO
import re
import os
path="C:/Users/nessr/Desktop/info projet"
os.chdir(path)
records=list(SeqIO.parse("plasmo.fasta","fasta"))
ATPM="[LIV]G[^P]G[^P][FYWMGSTNH][SGA][^PW][LIVCAT][^PD].[GSTACLIVMFY].{5,18}[LIVMFYWCSTAR][AIVP][LIVMFAGCKR]K"
patern=re.compile(r"(ATPM)")

def atp_binding(records,patern):

    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()

            iend=iend+50

            print(istart,m.group())

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

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

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

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

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

I would use something similar to this:

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

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

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


What methods should I use from PythonCyc API to query metabolites in BioCyc database? I am using PythonCyc API in order to write a query for metabolites in BioCyc. The purpose of this API is to communicate with the database software of BioCyc- Pathway Tools. Pathway Tools is in lisp therefore, PythonCyc creates a bridge between python and common lisp. To use this this API one must first create a PGDB object with a specified organism ID (orgid). In the example below I create a PGDB with orgid "meta". After this, I am able to call methods from PythonCyc with this object:

import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('sucrose')
print pathways

The query for the metabolite above, 'sucrose', provides a list a pathways:

[u'|PWY-7347|', u'|SUCSYN-PWY|', u'|PWY-7238|', u'|PWY-6527|', u'|PWY-5343|', u'|PWY-3801|', u'|PWY-7345|', u'|PWY-862|', u'|SUCUTIL-PWY|', u'|PWY66-373|', u'|PWY-621|', u'|PWY-822|', u'|SUCROSEUTIL2-PWY|', u'|PWY-6525|', u'|PWY-6524|', u'|PWY-5337|', u'|PWY-5384|']

However if I switch this metabolite name to a common amino acid name, such as 'valine':

import sys
import pythoncyc
#this creates PGDB object associated with meta(MetaCyc)
meta = pythoncyc.select_organism('meta')
#lists pathways of specified compound
pathways = meta.pathways_of_compound('valine')
print pathways

The query provides an error message that states "non coercible frame", which means it cannot find this ID. This is a common metabolite, just as sucrose, which should very well have pathway entries in database, however there are none found with this method. I have also tried synonyms of valine and other methods such as "all_pathways(cpds)" stated in API, which give me the same error message.
What is the reason for this? Does BioCyc not list any pathways for amino acids? Or am I using an incorrect method to access amino acid information?

Link to PythonCyc API

UPD: As suggested by Llopis, maybe the correct identifier is VAL?

I think there is no metabolite with such name in the database. I tried their search, the correct name seems to be L-valine. I expect that they use similar naming for other amino acids. If you cannot find something, you can use their compound search through the website.


Oxford nanopore promethion; does guppy3 separate fast5 into pass and fail? If not what does? just a quick question about fast5 pass fail folders. I have been working under the assumption that the fast5 pass fail folders I was given as raw data from our ONT vendor came from them doing base calling with guppy 3.something. Is that right? (The person that worked for the vendor has since moved on so I can't ask them - so my question is really does guppy3 divide fast5s into pass/fail?)
I further assumed that given I am now doing my own basecalling with guppy4.4 that I should use all the fast5s (both pass and fail) as some of the fails might become pass with the improved NN. Again, am I right? If not could someone please explain?
I'm going to also put this question on the ONT community site, FYI
Thanks,
Liam

Only if they've enabled the --qscore_filtering flag when using Guppy; this flag isn't enabled by default when running Guppy from the command line.
It's a good idea to use all fast5 reads for recalling, because accuracy improves over time, so some that previously failed may slide into the pass folder.
FYI, I'd be much happier if people asked their nanopore sequencing questions here. It's a publicly-accessible website which is indexed by search engines, so makes it easier to find answers, and easier to see the answers when found (i.e. no login needed).
Note that when recalling, you set the input and output directories (I think via the -i and -s flags). These should be different directories, and ideally different output directories for each recalling attempt, to avoid duplicate file entries. The original fast5 files shouldn't be modified by recalling as long as the input and output directories are different.


Spliced vs. unspliced ratios for transcripts in RNA-seq data Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?

I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)

Edit: (Additional information)

We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are.

Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions.

A more advanced approach I'm also considering is Salmon.

I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).

import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.

Once you've got these statistics you can calculate the "splicing index" as
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.


Why models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decrease I am reading Elowitz et al. (2002). When talking about intrinsic noise, i.e. noise due to microscopic events that govern which reactions occur and in what order, it is stated that:

  Models of stochastic gene expression predict that intrinsic noise should increase as the amount of transcript decreases.

Noise, defined in the article as the standard deviation divided by the mean, can be divided in intrinsic (due to microscopic events during transcription) and extrinsic (due to localization, state and concentration of other cellular components, e.g. polymerases).

I think the intrinsic noise should increase at low level of transcript amount because the average (mean) amount will be lower but deviation decrease is smaller than mean decrease. I don't know, however, if this is the actual reason (and why deviation decrease should be smaller than mean decrease) or if deviation decreases with the mean and there's another more likely explanation.

Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S. Stochastic gene expression in a single cell. Science 297, 1183–1186 (2002).

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

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

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


How to perform DE analysis for each sample I am new to R and biocondunctor. I have the normalized expression values for 20 samples for a disease and for 10 controls. I wanted to get the differential expressed values for each sample with all 10 controls.

How could I do that with limma.

Any help will be appreciated.

Edit: What I want to achieve is single-sample t-test between the control population and the single sample values to get DE values

Treat your disease samples as individual groups and then follow the normal routine in limma to use contrasts to compare two group (the control group versus individual disease samples). Note that you need to think long and hard about what these results then mean. Normally we compare groups because then the results should generalize to other samples. That will not be the case when you look at individual disease samples, so I suggest you do some clustering of the samples to see if you perhaps have any coherent clustering of the disease samples into subsets that could then be used as more useful groups in DE testing.


Getting hierarchy of cell populations with Drop-seq data I plan on collecting some drop-seq data from brain tissue (from mice and humans). Using only the 75-85% of genes that have 1:1 orthology between mice and humans, I'd like to cluster the cells by transcriptomic similarity and get a dendrogram to use for cell type classification. There are tons of tools for processing and clustering single cell RNA-seq data (https://www.scrna-tools.org/tools), but I get the impression that most of them were not designed with Drop-seq in mind. Anyone have any insight into the pros and cons of the best tools for Drop-seq analysis?

This article uses the freely available R package dropbead for filtering and then Seurat to perform a principal component analysis that groups together affine transcriptomes. It could be what you are looking for.


Variant calling without matched normal sample I have WGS .bam files for 3 patients (tumour and its matched derived model namely organoid) but I don't matched normal sample. If I call variants of each patients (tumour and its matched organoid), how I can use read counts at germline heterozygous positions estimated by GATK 3.2-2 to compensate for the absence of matched normal sample? I heard people use dbsnp VCF instead of the matched normal small variant VCF. But, I don't know start from where? which GATK function does that?

  I should mention, Calling SNV and indel in many tools returns vcf but
  I called copy number by varscan that did not return .vcf output so I
  am not sure what to do for CNV

Any helps please?

If you just want to filter out calls present in dbSNP then use:

java -jar GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R reference.fasta \
   -V patient.vcf \
   --discordance dbSNP.vcf \
   -o patient.filtered.vcf

--discordance will produce calls not present in dbSNP.


Accessing UCSC genome via ssh results in a validation error I am working on macOS High Sierra.

I am following the steps described here to lift an annotation over from one version of a genome to another. I am now just using the example genomes provided in the tutorial.

At this step:

time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
      -ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
  -bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
  -target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
  -query2Bit=../${query}/${query}.2bit \
  -querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1

I get the following error in the do.log file:

Warning: Permanently added 'localhost' (ECDSA) to the list of known hosts.
username@localhost: Permission denied (publickey,password,keyboard-interactive).
HgStepManager: executing from step 'align' through step 'cleanup'.
HgStepManager: executing step 'align' Thu Nov  8 14:27:32 2018.
Using localhost, /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.2bit and ../GCF_000004515.3_V1.1/GCF_000004515.3_V1.1.2bit
align: localhost does not have /Users/username/Documents/projects/main/output/software/liftover/data/genomes/GCA_000004515.3_Glycine_max_v2.0/GCA_000004515.3_Glycine_max_v2.0.ooc -- if that is not the correct location, please run again with -ooc.

I think the main problem is here:

username@localhost: Permission denied (publickey,password,keyboard-interactive).

This post describes this issue. However:

I do not have this file: ~/.ssh/authorized_keys. My files in ~/.ssh are: github_rsa github_rsa.pub known_hosts. One of the lines of known_hosts is localhost. Does this presumably allow me to ssh to my own laptop? I still don't understand why I need to ssh in this script if it's run locally.
PasswordAuthentication yes is already in /etc/ssh/sshd_config
If I run the command without the time command, I still get the same error. I am not prompted to give a password.

Could anyone help me overcome this problem, please?

You need to setup ssh keys.  All the commands in the script are going to use ssh to run the commands.  Note the tips here of how to setup ssh keys:

http://genomewiki.ucsc.edu/index.php/Parasol_job_control_system#SSH_keysssh key setup


Is there a way to assemble contigs starting from a specific sequence? My work involves searching for marker genes/fragments in metagenomic databases (like the Sequence Read Archive). Once I find these sequences, I would like to know more about the neighboring genomic region.

Is there a way I could assemble only sequences that create a contig which contains my region of interest? Contigs which don't contain this region are not useful to me. My organism of interest might represent a minority of the metagenome, and assembling everything in the dataset would use a lot of computing power.

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


Summing two columns where they share interval overlap I have two bed files where each column has a numeric score. I want to sum these numbers where there is overlap between the two bed files. How do I do that? Also, I am only interested in the score in specific genes.

Here is my setup:

import numpy as np
np.random.seed(0)
import pyranges as pr
import pandas as pd

entries = int(1e6)
cpg1 = pr.random(entries)
cpg2 = pr.random(entries)

cpg1.CPG = np.random.randint(100, size=entries)
cpg2.CPG = np.random.randint(100, size=entries)

regulatory_elements = pr.random(int(1e5), length=int(1e4))
regulatory_elements.Gene = np.arange(len(regulatory_elements), dtype=int)

So my data looks like:

cpg1
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | CPG       |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 8830650   | 8830750   | +            | 72        |
# | chr1         | 9564361   | 9564461   | +            | 57        |
# | chr1         | 44977425  | 44977525  | +            | 99        |
# | chr1         | 239741543 | 239741643 | +            | 54        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 29437476  | 29437576  | -            | 44        |
# | chrY         | 49995298  | 49995398  | -            | 43        |
# | chrY         | 50840129  | 50840229  | -            | 10        |
# | chrY         | 38069647  | 38069747  | -            | 89        |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 1,000,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

regulatory_elements
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Gene      |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 217921634 | 217931634 | +            | 0         |
# | chr1         | 166226804 | 166236804 | +            | 1         |
# | chr1         | 170688210 | 170698210 | +            | 2         |
# | chr1         | 57958563  | 57968563  | +            | 3         |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 35870043  | 35880043  | -            | 99996     |
# | chrY         | 35057634  | 35067634  | -            | 99997     |
# | chrY         | 52718975  | 52728975  | -            | 99998     |
# | chrY         | 19196507  | 19206507  | -            | 99999     |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 100,000 rows and 5 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

(This is a q from privateish correspondence).

If you want to sum two columns based on overlaps, converting them to RLEs is a good idea.

# get the cpg scores only where there is overlap with the genes you are interested in
cr1 = cpg1.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")
cr2 = cpg2.join(regulatory_elements).new_position("intersection").drop(like="_(a|b)")

# convert these to run length encodings, using the score from the CPG column
rle1 = cr1.to_rle("CPG", rpm=False)
rle2 = cr2.to_rle("CPG", rpm=False)

rle_sum = rle1 + rle2

# convert the rles back to ranges
# (I will give the to_ranges function an option dtype eventually)
range_sum = rle_sum.to_ranges().apply(lambda df: df.astype({"Start": np.int32, "End": np.int32}))

# add the gene names to the scores
gr = range_sum.join(regulatory_elements) # .drop(like="_b")

Now the result is:

print(gr)
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 15623     | 15723     | 77.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 16034     | 16134     | 84.0        | +          | 9886      | 19886     | +            | 3714      |
# | chr1         | 36588     | 36688     | 2.0         | +          | 36479     | 46479     | +            | 2871      |
# | chr1         | 37271     | 37371     | 30.0        | +          | 36479     | 46479     | +            | 2871      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59312851  | 59312951  | 28.0        | -          | 59308107  | 59318107  | +            | 98599     |
# | chrY         | 59332843  | 59332943  | 82.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59337592  | 59337692  | 41.0        | -          | 59332664  | 59342664  | -            | 99215     |
# | chrY         | 59341719  | 59341819  | 15.0        | -          | 59332664  | 59342664  | -            | 99215     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 668,317 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

print(gr[gr.Score > 100])
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Score       | Strand     | Start_b   | End_b     | Strand_b     | Gene      |
# | (object)     | (int32)   | (int32)   | (float64)   | (object)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------|
# | chr1         | 96874     | 96974     | 116.0       | +          | 87668     | 97668     | -            | 7391      |
# | chr1         | 96874     | 96974     | 116.0       | +          | 89794     | 99794     | -            | 6357      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 744041    | 754041    | +            | 1150      |
# | chr1         | 749686    | 749697    | 180.0       | +          | 749277    | 759277    | -            | 6395      |
# | ...          | ...       | ...       | ...         | ...        | ...       | ...       | ...          | ...       |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59052942  | 59062942  | -            | 99957     |
# | chrY         | 59060667  | 59060767  | 196.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59063570  | 59063576  | 105.0       | -          | 59059904  | 59069904  | +            | 98401     |
# | chrY         | 59269010  | 59269085  | 152.0       | -          | 59264704  | 59274704  | -            | 99499     |
# +--------------+-----------+-----------+-------------+------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 106,643 rows and 9 columns from 25 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.


Index VCF file for rapid access by rsID? With tabix one can index a VCF file for rapid variant retrieval based on genomic position. I'm wondering if there are any tools that will index a VCF file to allow rapid retrieval using rsIDs and/or other metadata? I'm aware of awk/grep/vcftools one-liners for this purpose, but I'd like to avoid scanning a huge VCF each time I need to retrieve the coordinates of a new batch of rsIDs.

I am not aware of any widespread index formats for query-by-QNAME for BAM or query-by-ID for VCF. I thought there were historical samtools{-devel,-help} threads to the effect of “in principle one could index a named-sorted BAM file, but there's insufficient demand so no-one's ever implemented it”, but I can't find any just now.

Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.

If you are using dbSNP, probably you have a dbSNP database or API that you can use to query by rsID and get genomic position back. So you can retrieve your VCF contents by rsID in two stages: do that dbSNP query, and use its results to find the variants in your VCF file by position.


Duplicate gene symbol handling in GEO gene expression data I have downloaded a gene expression data from GEO database (GSE3268) in which in some of its rows there are duplicate gene symbols.
For example TP53 exists in two rows with different expression values and different GenBank Accession IDs.
They have the following Nucleotide Title:

Human p53 cellular tumor antigen mRNA, complete cds
Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA

I want to map these data with edge list of gene symbols to make a network. Which expression value should I consider for such duplicate gene symbols for mapping?

Since this is a microarray (U133A), what you're seeing is that each probe is associated with one or more transcripts. The general strategy when dealing with microarrays is to use RMA (or fRMA) to summarize to probesets and then collapseRows to summarize to gene-level information. Since this is all typically done in R, you can find some further discussion on the bioconductor site. I would suggest going through this process and then only using the result of that in building your network (whether you're happy with sticking to transcript-level quantifications or gene-level quantifications for building your network is up to you).


Simplest way to work out structural variant type? In VCF 4.2, a structural variant (SV) can be described with the BND keyword in SVTYPE. For example, the following example is an insertion (from https://samtools.github.io/hts-specs/VCFv4.2.pdf):

#CHROM    POS ID REF ALT                      QUAL FILTER INFO
2 321682  bnd V  T   ]13 : 123456]AGTNNNNNCAT 6    PASS   SVTYPE=BND;MATEID=bnd U
13 123456 bnd U  C   CAGTNNNNNCA[2 : 321682[  6    PASS   SVTYPE=BND;MATEID=bnd V

The keyword can be used to describe any SV event.

I would like to work out how to find out if an event is an insertion, deletion, duplication, inversion, etc. The algorithm/program will tell me the example above is an insertion.

I'm not able to find a tool that can do the conversion for me.

PS: The program I'm using is GRIDSS. This program would only generate BND in the SVTYPE field. It's not strictly a bug (BND can be used for any SV arrangement), but it's annoying in analysis.

GRIDSS's author had posted his own solution on Github. The code is:

simpleEventType <- function(gr) {
    return(ifelse(seqnames(gr) != seqnames(partner(gr)), "ITX", # inter-chromosomosal
           ifelse(gr$insLen >= abs(gr$svLen) * 0.7, "INS",
           ifelse(strand(gr) == strand(partner(gr)), "INV",
           ifelse(xor(start(gr) < start(partner(gr)), strand(gr) == "-"), "DEL", "DUP")))))
}


Finding transposable elements using RepeatMasker I'm using RepeatMasker to detect, classify the Transposable elements. My Input is a eukaryotic non-reference genome.
I made a run via RepeatMasker many times to Mask the TEs, but return 0 Annotation tables.
Further, I used a different -species option each time, but Annotation results contain only short and simple repeats.

Do I have to create a unique library to detect them?

The command I used:

./RepeatMasker -no_is -noint - species mammal Genome.fna

Results:

sequences:         32573
total length: 2004063690 bp  (1981588036 bp excl N/X-runs)
GC level:         41.28 %
bases masked:   33282890 bp ( 1.66 %)
==================================================
               number of      length   percentage
               elements*    occupied  of sequence
--------------------------------------------------
SINEs:                0            0 bp    0.00 %
      ALUs            0            0 bp    0.00 %
      MIRs            0            0 bp    0.00 %

LINEs:                0            0 bp    0.00 %
      LINE1           0            0 bp    0.00 %
      LINE2           0            0 bp    0.00 %
      L3/CR1          0            0 bp    0.00 %

LTR elements:         0            0 bp    0.00 %
      ERVL            0            0 bp    0.00 %
      ERVL-MaLRs      0            0 bp    0.00 %
      ERV_classI      0            0 bp    0.00 %
      ERV_classII     0            0 bp    0.00 %

DNA elements:         0            0 bp    0.00 %
     hAT-Charlie      0            0 bp    0.00 %
     TcMar-Tigger     0            0 bp    0.00 %

Unclassified:         0            0 bp    0.00 %

Total interspersed repeats:        0 bp    0.00 %

Small RNA:            0            0 bp    0.00 %

Satellites:           1          267 bp    0.00 %
Simple repeats:  663679     27953718 bp    1.39 %
Low complexity:  123705      6154714 bp    0.31 %
==================================================

You need to model the repeats in your de novo genome.

See:
http://www.repeatmasker.org/RepeatModeler/ and
https://www.biostars.org/p/154290/ (from previous answer)

You first build a database named "name_of_your_database" (insert whatever you want to call it) and run it on the fasta file of your genome (e.g. your_genome.fasta).

Make you find the right paths for the software.

## build the database
BuildDatabase -engine ncbi -name name_of_your_database your_genome.fasta

## run the modeler (in this case using ncbi)
RepeatModeler -database name_of_your_database -engine ncbi


Fast filtering of intervals not falling within a certain distance from known genes I would like to filter a bed file with intervals, ie. in the format of:

chr1    13800   14301
chr1    15500   16001
chr1    19400   19901
chr1    22800   23301

In particular, I want to filter out intervals that fall far from known genes, let's say > 1KB.

I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from Homo_sapiens.GRCh38.93.gtf by filtering for gene coordinates:

1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";

However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

Convert GTF annotations to a sorted BED file of genes:

$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

Sort your intervals, if unsorted:

$ sort-bed intervals.unsorted.bed > intervals.bed

Filter intervals against 1kb-padded genes:

$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed

We first --merge genes to simplify dealing with overlapping genes and further speed downstream filtering.

The file intervals.filtered.bed will contain intervals from intervals.bed which lie outside of a 1kb-window around merged gene space.

To run all of this as a one-liner (assuming use of the bash shell):

$ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed

This creates the same output result, but avoids creation of intermediate files and so should run even faster.

Note: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using --range takes care of ranges that would have a negative start position, setting their start coordinates to 0. Bounds aren't really relevant here, in any case, because the bedops --not-element-of step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.


What are the ways to keep track of branches in the analysis? I'm going through an RNA-seq pipeline in R/Bioconductor and want to try multiple parameters at subsequent steps, for example, running clustering with different settings, running RegressOut or not on unwanted effects etc. That's a lot of "versions", even if I don't do combinations of these steps.

How can I keep track of this, and my conclusions?
Not necessarily want to save the results.

Save the different scripts with git (seems overkill)
Make notes in the script itself

Save the different scripts with git (seems overkill)

Whoa. I did an actual double take when reading this:1 it’s the opposite of overkill.

Version controlling your scripts (using Git or something similar) is the absolute minimum, and should become completely automatic. For every new project I begin, one of the very first steps is to issue the git init command, and to set up a remote repository (on Github).

To keep track of different analyses, I use a combination of the following approaches:

Write reusable functions/scripts and parametrise. The parameters are kept either inside the script itself (that then calls the relevant function repeatedly), or in a Makefile (I recommend Snakemake).
Document the alternative analysis approaches; once again, this could be a Makefile with different rules for alternative analyses, or a set of notebooks (via R Markdown).
Have different Git branches for mutually exclusive approaches. At the end of the analysis one of these branches gets merged into master, and published. If I want to publish several analysis approaches, I merge all these branches into master, and use approaches (1) or (2) enable them simultaneously.

In fact, I recommend creating a Makefile for every analysis; I have found that this is the most practical, self-documenting way to run an analysis. It most closely resembles a wet-lab lab notebook. The advantage over a single R Markdown document is that rerunning just parts of the analysis can be completely automated, and dependencies in the workflow are apparent from the dependencies of the Makefile rules. This is much harder in R Markdown.

Some time ago I create an example analysis workflow to show how this can be structured. Nowadays I would use Snakemake instead of GNU make.

Regarding your other point:

  Make notes in the script itself

“Notes” are a dangerous beasts: documentation is important, but experience shows that it’s sometimes very hard to keep documentation synchronised with the code. There is no mechanism that ensures that documentation and code actually agree. Differences between presumed analysis and actually executed analysis can become very problematic.

People therefore prefer using self-documenting code as much as possible; that is: writing code so that its meaning becomes immediately clear, without comments, even to somebody who hasn’t worked on the code before. Doing this well is hard and takes practice, but improves overall code quality. Once again, using a Makefile helps here because the dependencies between the rules are self-documenting the kind of analysis that was performed.

Jeff Atwood has written two seminal essays on this subject:

Coding without comments, and
Code tells you how, comments tell you why.

They are two of the best pieces of advice on programming that I could give.

1 To emphasise: take a look at the edit history of this answer.


WGBS vs ONT, which one should I trust? First of all let me apologize if this is not the place for this question. While it can be quite broad, I will try to make it as specific as I can.
Context
I am testing the correlation between Whole Genome Bisulfite Sequencing and Oxford Nanopore methylation calls for one individual (mammal, blood tissue) with replicates at all CpG sites. In the (very) best cases I find a close to 0.9 correlation. Yet, for the part that does not correlate I would like to know which technology is closer to the "biological truth".
My strategy
My plan would be to first find CpG sites/short CpG regions showing differential methylation states between WGBS and ONT. Then I would validate the methylation states of these sites using a third party technology such as pyromark.
What I am asking to the SE community

First, has anyone done this before ? I would love to hear about your experience

Do you see some critical aspects I should be watching out for ?

Last if you have comments about the strategy I propose here, you are welcome

Methylation can be cell-specific, which makes it difficult to evaluate accuracy on a bulk-cell level (even within the same tissue). How can you tell that the differences you're seeing are due to platform differences, or due to biological variation?
I find that adding more haystacks doesn't help much in working out the truth of a dataset. If you want to investigate biological truth, it would be better to create a biological system with a known methylation (or demethylation) pattern and test that.


How to concatenate two fasta sequences by py I have three backbone and i want to concatenate 70 sequence into these backbone.
such like:
fasta file 1:
Backbone1: ACTGCGTATGTCAGCATATAGACCGCCT
Backbone2: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTC
Backbone3: GTCATGCTACGTGCTACGTGCTCGTCGCGC

fasta file 2:
sequence_1:CGCTAGTGTAGCTGCTACTCCG

output sequence:
Backbone1+sequence_1: ACTGCGTATGTCAGCATATAGACCGCCTCGCTAGTGTAGCTGCTACTCCG
Backbone2+sequence_1: CGTAGCTAGCTAGCTCGTGTCGTGATGCTACTCCGCTAGTGTAGCTGCTACTCCG
Backbone3+sequence_1: GTCATGCTACGTGCTACGTGCTCGTCGCGCCGCTAGTGTAGCTGCTACTCCGs

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

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

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


At what processing step should library strandedness type be taken into account? Suppose I have single-end RNA-seq data for which the reads in the fastq file are reversed with respect to the original extracted RNAs.

Suppose I have the following workflow:

Map the reads on the reference genome, sort and index the resulting sam output.
Use the bam file to make a bigwig file showing reads mapping on the plus and minus strands of the chromosomes.
Use the bam file and an annotation file to quantify the gene expression.

Steps 2 and 3 depend on step 1, but are not dependant on one another.

I'm sure that strandedness needs to be taken into account no later than at step 3.

The questions are the following:

Is it common to have mapping tools with an option to "correct" for strandedness?
Assuming the bam file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Regarding this second point, I'm split between two attitudes:

There should be no correction, the bigwig files are there to represent the "raw" mapping results.
The bigwig files are there to give an idea of the origin and abundance of the originally extracted RNAs, so a strandedness correction should be used if necessary.

Which one is the most common attitude?

If the second one is preferred, then wouldn't it be "safer" / "cleaner" to handle strandedness as early as possible? Possibly at the mapping, or even using a pre-processing step?

I usually take into account library strandness during the mapping to the reference step (but I never worked with bigwig files myself).

Assuming RNA-seq data from Illumina, you can use Hisat2 for alignment to the reference genome, and you have the option --rna-strandness to specify the strand-specificity information (default unstranded). This is particularly important when you wish to have XS strand attribute in your bam file for easier compatibility with downstream software (especially important in my experience if you are planning to use some utility from Cufflinks for your step 3).


Genome annotation of a Bacillus strain I finished the sequencing of the genome of a biocontrol activity strain Bacillus amyloliquefaciens by Illumina hiseq. I have already gone through some papers but there were no details of processing in it.

Is anywhere I can find some resources to get my research going on to genome annotation, compare it with other microbial genome data, and analyze non-ribosomal peptide synthetase, find out more genomic information about biocontrol activity?

For example, how to locate 16S DNA and some lipopeptide antibiotic genes (fungycin gene, surfactin gene) in my bacterial genome, and compare it with the reference strain and draw a phylogenetic tree. And how to use my bacterial genome and 2 reference genomes to draw a genomic circle diagram.

You said that you want to compare with other microbial genome data. Have you already annotated your genome? If not, you might consider a pipeline I've developed, I can gladly help if you face some problems using it.

It is called bacannot, it does a complete annotation using Prokka, Kofam, CARD, VFDB, and other resources. It might help you to have a first scratch of your data since it also collapses annotation information in a genome browser


Human Cell Atlas - Data availability A news item from 2017-10-18 on the website of the Human Cell Atlas states:

  In addition, the consortium today also announced the impending release of gene expression profiles from the first one million immune cells collected under the HCA, toward an initial milestone of collecting at least 30 million cells representing several tissues and organs for the atlas’ first draft. These data, to be posted on an online repository by early November, will be freely available for researchers’ use. (https://www.humancellatlas.org/news/14)

As of now, Feb 2018, I cannot find these data anywhere.
Does anybody know where I can download the gene expression data of these 1 million immune cells or when it will be eventually published?

There is now a post on Biostars:

  We are starting to define the release process for the Human Cell
  Atlas.

  The first couple of HCA releases will be primarily single-cell RNA-Seq
  and we are interested to know what you would expect our releases to
  contain

  Please fill out our survey if you would like to help us define our
  releases

  https://www.surveymonkey.com/r/S6Y753P

And I suppose that answers the original question: the data has not been released yet and probably won't be for a while considering they are still discussing how to make it available.

Update (from the same post):

  The first production release of the DCP is happening in June.

  The first data release of HCA is targetting Q4 2018.


Does UCSC definition for telomere making sense? UCSC provides a database of telomere regions for each chromosome:

  ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/gap.txt.gz

It reported the telomere section in chr21 as:

grep chr21 gap.txt | grep telomere
585 chr21   0   10000   1   N   10000   telomere    no
952 chr21   48119895    48129895    515 N   10000   telomere    no

Telomere is the first 10000 bases in chr21?? But ... the region are all "N" in the reference genome!

How does UCSC confidently give out telomere regions if they don't even know anything about the sequences?

It seems so, yes.

UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.

Taken from Wikipedia:

A telomere is a region of repetitive nucleotide sequences at each end of a chromosome, which protects the end of the chromosome from deterioration or from fusion with neighboring chromosomes. For vertebrates, the sequence of nucleotides in telomeres is AGGGTT, with the complementary DNA strand being TCCCAA, with a single-stranded TTAGGG overhang. This sequence of TTAGGG is repeated approximately 2,500 times in humans. In humans, average telomere length declines from about 11 kilobases at birth[ to less than 4 kilobases in old age, with the average rate of decline being greater in men than in women.

2,500 x 6 (TCCCAA) = 15,000
(approx. 10K)

Two interesting additional facts:

The minimum telomere length needed to ensure human telomere protective stability in white blood cells is 3.81 kb
Mouse has a 100Kb telomere representation in USCS


Error in creating a volcano plot in MATLAB I am a complete newbie to MATLAB and the first task I have is to create a volcano plot. I have been using this link to understand about it and get started.

I tried to run it on dummy values -

a=[1 2 3]

b=[4.6 2.7 4.5]

c=[0.05 0.33 0.45]

And then I ran -

SigStructure = mavolcanoplot(a, b, c)

My understanding is that a represents the gene expression values for condition 1, b for condition 2, and c is the list of p-values for the 3 data points in a and b.

However running this code gives me the error -

Index exceeds matrix dimensions.

Error in mavolcanoplot (line 127)
appdata.effect = X(paramStruct.goodVals) - Y(paramStruct.goodVals);

Error in volc (line 4)
SigStructure = mavolcanoplot(a, b, c)

Can anyone explain where I am going wrong?

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

Update: here is the actual script:

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

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

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


Rapid metagenomics classifiers on long read data I recently used the minION (Nanopore, 9.4 flow cell, RAD001 kit) to generate a metagenome out of environmental samples.
Passed reads weren't brilliant (196, average 1,594bp lenght), but working with centrifuge the classification outputs turned out to have quite low hitLength to queryLength scores (average 2%, max 14%). This plus score values don't give me a lot of confidence towards the results I got.

Has anyone else used centrifuge and experienced the same?

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

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

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


Is there a Python program that is able to show convergent and divergent evolution on a phylogenetic tree? I was wondering if there was a Python program that detects convergent and divergent evolution on a tree. I am also curious if it is just better to look at a phylogenetic tree and conclude what is convergent or divergent evolution than what a program tells you.

Well in the absence of a definition of 'divergent evolution' I would simply use a definition of allopatry. If the divergence is congruent with geographic distance or geographic foci this provides some proof of divergent evolution.
In this context you are looking for something like the biogeographic BAyesian analysis such as BayArea or geographic add-ons to Beast
There isn't a singular Python library that will handle this but there is some functionality with Biopython here
from Bio import Phylo


How to find GEO data sets using Drug Bank ID and bioDBnet I want to find all experiments in GEO that are associated with a drug (for example tolvaptan). Is there any quick and scalable way to to this? I want to query more than 100 drugs.

I tried to use bioDBnet to map Drug Bank ID to look up data sets, but I don't know what output format should I select.

You could take a look at the rentrez package for R. You can use the entrez_search function to get the GEO IDs for the studies related to your query, and then the entrez_summary function for the results. So, for example:

x <- rentrez::entrez_search(db="gds", term=paste0("tolvaptan AND gse[ETYP]"),
                            use_history = T)
res <- rentrez::entrez_summary(db="gds", id=x$ids)
str(res)

Where res will be a list containing objects of the esummary class. Then, you can use the following line to get the GSE accession numbers related to tolvaptan:

lapply(res, `[[`, "gse")

So you could put this function in a loop replacing "tolvaptan" with your drugs:

drugs <- c("drugA", "drugB", "drugC")
for(drug in drugs){
    x <- rentrez::entrez_search(db="gds", term=paste0(drug, " AND gse[ETYP]"))
    res <- rentrez::entrez_summary(db="gds", id=x$ids)
    lapply(res, `[[`, "gse")
}

Take a look at the rentrez documentation for more details.


Splitting data table up based on separate file of gene names I am very new so bare with me. I have a file called data (.csv that I turned into a data table) with 2 columns (geneid, normalized counts) and about 9000 rows. I have a second file called Xgenes of genes that I am interested in pulling out of the first table (1 column- geneid, about 50 rows).
Is there a way to make 2 new files- I want one with all the geneid/normalized counts minus the ones in the second file (so minus about 50 genes). Then I want to make another file that only has the geneid/normalized counts of the 50 genes in the second file.
Currently I am doing this manually which does indeed work but if there is a way to do it through R it would make my future analysis much easier.
Thanks!

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

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

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

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


How to convert DNAbin to FASTA in R? I am trying to convert DNAbin files to fasta format. Rationale: I want to use the fasta file to calculate non-synonymous/synonymous mutation rate.

my_dnabin1 is a DNAbin file of 55 samples and I am using the following code to convert it into a fasta file.

dnabin_to_fasta <- lapply(my_dnabin1, function(x) as.character(x[1:length(x)]))
This generates a list of 55 samples which looks like:

$SS.11.01

[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

$SS.11.02

[1] "t" "t" "a" "c" "c" "t" "a" "a" "a" "a" "a" "g" "c" "c" "g" "c" "t" "t" "c" "c" "c" "t" "c" "c" "a" "a"
 [27] "c" "c" "c" "t" "a" "g" "a" "a" "g" "c" "a" "a" "a" "c" "c" "t" "t" "t" "c" "a" "a" "c" "c" "c" "c" "a"

and so on...

However, I want a fasta formatted file as the output that may look something like:

>SS.11.01 ttacctga

>SS.11.02 ttacctga

The following should work:

cat(file="out.fa", paste(paste0(">",names(dnabin_to_fasta)),
          sapply(dnabin_to_fasta, paste, collapse=""), sep="\n"), sep="\n");

A bit of explanation, starting from the middle outwards:

Collapse the letter arrays inside dnabin_to_fasta into single strings
Append ">" to the front of the names of each sequence
Concatenate the results of steps 1 & 2 together, separating with a line break
Write each concatenated string to a file, separating by a line break


Is there a way to get the code from "github.mit.edu"? everyone.
I am a postdoc who just begins to analyze single-cell RNA-seq data.
Nowadays, I found a really interesting paper (https://doi.org/10.1016/j.celrep.2018.10.047, PMID: 30404002). So, I decided to study bioinformatics by drawing the same graphs presented in the paper during the corona outbreak season. However, it was turned out that I can't get any access to the website in which code is deposited. Here is the address for the website.

Do you guys know how to get the code to learn bioinformatics and make graphs?

Thank you.

It is very strange that they put their code somewhere cannot be easily accessed. If they just put it under a normal Github repo, there would not be such a problem. One way to get the code is simply to email the corresponding author. People are usually responsive if you are interested in their work.


bed file to gtf conversion using awk I have a bed file like this:

chr1    12227   12612   LOC102725121    0   +
chr1    12721   13220   LOC102725121    0   +
chr1    12227   12612   DDX11L1 0   +
chr1    12721   13220   DDX11L1 0   +
chr1    14829   14969   WASH7P  0   -
chr1    15038   15795   WASH7P  0   -
chr1    15947   16606   WASH7P  0   -
chr1    16765   16857   WASH7P  0   -
chr1    17055   17232   WASH7P  0   -
chr1    17368   17605   WASH7P  0   -

I am trying to convert into gtf file for downstream analysis. Gtf file should look like this:

chr1   myIntrons   feature 3073253   3074322   .   +   .     gene_id "intron1"; transcript_id "intron1"
chr1   myIntrons   feature 3075858   3079636   .   +   .     gene_id "intron2"; transcript_id "intron2"

For the gene_id part in the 9th column in the double quotes I want my 4th column and in transcript_id I want combination of column1,column2,and column3 so that my output looks like this:
Showing for the first line

chr1    myIntrons   intron  12228   12612   .   +   .   gene_id "LOC102725121"; transcript_id "chr1:12227-12612"

I am using this code but under double quotes it is not expanding $1 $2 etc

awk -v OFS="\t" '{print $1,"myIntrons","intron", $2+1, $3, ".", $6, ".", "gene_id $4; transcript_id $1:$2-$3"}' file.bed| head

Instead of the "gene_id ..." bit, you want sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"", $4, $1, $2, $3).


Cutting the sequence into several sequences with the information of a dataframe I have a fasta file with several 120-concatenation protein sequences. I also have a data frame with 120 names of proteins in column 1 and their length in column 2. Using this dataframe I want to create a new fasta file per sequence in which I get each protein independent with the name of the protein and the bacteria.

Here my fasta file.

Bacteria_1
Seq('GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...IEE', SingleLetterAlphabet())

Bacteria_2
Seq('GRRKDSVARVYLRPGKGEIVVNKRPVEDYFGRETLKMVLRQPLELTESLDKFDI...YK-', SingleLetterAlphabet())

Here my dataframe

      Names  LEN
0     Ribosomal_S9:  121
1     Ribosomal_S8:  129
2    Ribosomal_L10:  100
3             GrpE:  166
4           DUF150:  141
..              ...  ...
115      TIGR03632:  117
116      TIGR03654:  175
117      TIGR03723:  314
118      TIGR03725:  212
119      TIGR03953:  188

[120 rows x 2 columns]

The name of the genes and their length is the same for all my bacteria so this df is applicable for all sequences in my fasta file.

I am trying to get new fasta files like this:

>Bacteria_1 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_1 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_1 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

The second fasta file

>Bacteria_2 Ribosomal_S9:
GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDV...

>Bacteria_2 Ribosomal_S8:
KDYKDFFPVPTLAAVARAYMTAGKGNLGQYDVGRRKAAVARAYMTAGKGN...
...

>Bacteria_2 TIGR03953:
AJEFISUBEVFCVALSKDMCFPADNSCF0AISCIABNSCNAISCN...

We wouldn't write the code for you from scratch (thats your job), but its an easy problem because in Python any string is an "array" or "dictionary". Thus

bacteria1 = 'GRRKAAVARAYMTAGKGNVTVNKKDYKDFFPVPTLQYKITQPFAITETLGQYDVIEEXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX'
print (bacteria1[0:120])

This will define the first gene in your list, remembering the array starts at 0 whilst your gene starts at position 1 (so you need to subtract 1). You simply parse the input list so the gene lengths form a dictionary and loop through the dictionary of gene lengths ensuring the last length is added to the current one and then issuing a write command.

You don't really need biopython in this instance, unless it will parse the gene length file you need to open (I'm not aware it can).

BTW the format you supplied isn't fasta it looks like a Biopython statement, again in this problem I wouldn't bother because Python is perfectly capable of splicing a string.


5'UTR and 3'UTR annotation in yeast I am working on a project in which I need to compute several parameters (such GC content and length) of 5'UTR and 3'UTR sequences of Saccharomyces cerevisiae yeast genes.

The problem is finding a proper annotation for these regions in yeast. I have tried with BiomaRt and these sequences are not available for yeast. I have also tried with UCSC table browser to obtain a BED file with 5'UTR or 3'UTR sequences but it is returning only coordinates of tRNAs, pseudogenes, ncRNAs.

In summary, is there a straight-forward way to obtain an annotation for 5'UTR and 3'UTR regions in Saccharomyces cerevisiae yeast?

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

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

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

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

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

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


How can I remove (non-trivial) duplicates from a VCF file? This is related to the question I asked here. Consider a vcf file that contains duplicate variants, but where the duplicates aren't simply the same thing in the same notation but instead one is a subset of the other. For example:

##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

These the two variants are actually the same. They result in exactly the same genotype. Changing  AACAC to AATAC at position 529514 just means change C to T at position 529516.

Is there any tool that can detect such duplicates and remove them? I tried vcfuniq from vcflib, but that doesn't seem to recognize this as a duplicate. I think it only looks at the 1st 4 fields and only considers duplicates those variants with exactly the same values in the 1st 4 fields:

$ ./bin/vcfuniq test.vcf
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  .   AACAC   AATAC   .   PASS    .   GT  0/1
chr12   529516  .   C   T   .   PASS    .   GT  0/1

However, as explained in the linked question, EBI's vcf_validator considers this invalid. And it doesn't really make sense to have these duplicates in any case, so is there any way I can detect and remove them? Preferably an existing tool, but I am open to scripting solutions as well.

This is further complicated by cases like this one:

##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC       0.00    .       .     GT    0/1
chr12   529516  529516  C       T       .       PASS    .     GT    0/1

Unfortunately, this one won't be caught by the approach in Daniel's clever script:

$ cat test2.vcf | foo.py
##fileformat=VCFv4.1
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529514  529514  AACAC   AAT,AATAC   0.00    .   .   GT  0/1
chr12   529516  529516  C   T   .   PASS    .   GT  0/1

It turns out that bcftools can do this (tested on bcftools-1.8), if you give it the reference genome to test against:

$ bcftools norm -d none -f hg19.fa test.vcf
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -d none -f hg19.fa test.vcf; Date=Wed Feb 27 16:08:44 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  .       C       T       .       PASS    .       GT      0/1
Lines   total/split/realigned/skipped:  2/0/1/0

For the more complex case of the multi-allelic variant in the second VCF example from the question, you can run it through bcftools twice. Once using norm to left-align and split multi-allelic variants, and then again to remove the duplicates:

$ bcftools norm -m -any -NO z -O v -o - ~/test2.vcf |
    bcftools norm -d none -f hg19.fa
Lines   total/split/realigned/skipped:  2/1/0/0
##fileformat=VCFv4.1
##FILTER=<ID=PASS,Description="All filters passed">
##reference=foo
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr12>
##bcftools_normVersion=1.8+htslib-1.8
##bcftools_normCommand=norm -m -any -NO z -O v -o - test2.vcf; Date=Wed Feb 27 18:18:32 2019
##bcftools_normCommand=norm -d none -f hg19.fa -; Date=Wed Feb 27 18:18:32 2019
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  Sample1
chr12   529516  529514  CAC     T       0       .       .       GT      0/1
chr12   529516  529514  C       T       0       .       .       GT      0/0
Lines   total/split/realigned/skipped:  3/0/2/0


How to create a list of differentially expressed (DE) genes after normalization with RUVSeq? I am using edgeR to perform differential expression (DE) analysis on a set of RNA-seq data samples (2 controls; 8 treatments). To correct for batch effects, I am using RUVSeq.
I am able to get a list of DE genes without normalization:
x <- as.factor(rep(c("Ctl","Inf"),c(2,8)))
set <- newSeqExpressionSet(as.matrix(counttable),phenoData=data.frame(x,row.names=colnames(counttable)))
design <- model.matrix(~x, data=pData(set))
y <- DGEList(counts=counts(set), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
top <- topTags(lrt, n=nrow(set))$table
write.table(top, paste(OUT, "DE_genelist.txt", sep=""))

Then immediately after creating the "top" object, I use RUVg to normalize:
# [...]
top <- topTags(lrt, n=nrow(set))$table
empirical <- rownames(set)[which(!(rownames(set) %in% rownames(top)[1:5000]))]
ruvg <- RUVg(set, empirical, k=1)
write.table(ruvg, paste(OUT, "DE_RUVg_genelist.txt", sep=""))

And I get the error:
Error in as.data.frame.default(x[[i]], optional = TRUE) :
  cannot coerce class ‘structure("SeqExpressionSet", package = "EDASeq")’ to a data.frame

I am not sure how to print the list of normalized results like I can with the unnormalized data. Ideally, I would get a file with the same format as the edgeR output (as a .csv or .txt file):
"logFC" "logCPM" "LR" "PValue" "FDR"
"COBLL1" -2.150 4.427061248733 75.0739519350016 4.53408921348828e-18 9.51203608115384e-15
"UBE2D1" -2.178 3.577168782408 74.9346752854903 4.86549160161322e-18 9.51203608115384e-15
"NEK7" -2.404 4.020072739285 72.6539117671717 1.54500340443843e-17 2.71843349010941e-14
"SMC6" -2.300 5.674738981329 61.8130019860261 3.7767230643666e-15 3.4974443325016e-12

How can I get a list of genes as an output after normalization with RUVSeq?

You do the normalization before running your edgeR. The purpose of RUVg is to remove "Remove Unwanted Variation Using Control Genes". In your code, you ran edgeR and then normalize the data using RUVg, which is only going to return you the normalized counts.
Using the example dataset in vignette:
library(RUVSeq)
library(zebrafishRNASeq)
data(zfGenes)
filter <- apply(zfGenes, 1, function(x) length(x[x>5])>=2)
filtered <- zfGenes[filter,]
genes <- rownames(filtered)[grep("^ENS", rownames(filtered))]
spikes <- rownames(filtered)[grep("^ERCC", rownames(filtered))]

x <- as.factor(rep(c("Ctl", "Trt"), each=3))
set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(x, row.names=colnames(filtered)))
set <- betweenLaneNormalization(set, which="upper")

set1 <- RUVg(set, spikes, k=1)

You can look at it, it's an expression set with counts etc, not results:
set1
SeqExpressionSet (storageMode: lockedEnvironment)
assayData: 20865 features, 6 samples
  element names: counts, normalizedCounts, offset
protocolData: none
phenoData
  sampleNames: Ctl1 Ctl3 ... Trt13 (6 total)
  varLabels: x W_1
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation:

You run edgeR now on the results of RUVg:
design <- model.matrix(~x + W_1, data=pData(set1))
y <- DGEList(counts=counts(set1), group=x)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
fit <- glmFit(y, design)
lrt <- glmLRT(fit, coef=2)
topTags(lrt)


How to score how densely bases are distributed along a given oligo? Suppose I have a short chain string of oligos, and I want to rate them from 0 to 1 based on how clustered the "A"s are in respect to the chain. For example:

"AAAAACCCCC" #should score as 1 for being as clustered as it possibly could
"ACACACACAC" #should be scored as 0 for being as thinly distributed as it could

Current idea is to strip both ends of A's, then remove first and last characters, then

gaps = ''.join(str(x) for x in deck).split("A")
listofCs = [x for x in gaps if x]

to join it all together as a list of how big the gaps are (i.e. the "C"s), and simply count the average size

average = sum(map(len,listofCs))/len(listofCs)

But I feel like there should be a much more elegant solution to this. It feels messy. How can I improve the scoring?

You could count the number of transitions between an A and a 'not A', and divide this number by the total number of As.  You might also need to scale by the sequence length, if you have different length sequences that you want to compare.

Edit: also, there might be a correction to make if one or both ends of the sequence are an A


How can I colour boxes in Gviz AnnotationTrack in R? I'm learning the Gviz bioconductor package, I generate a plot as follows:

library(Gviz)
track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'))
pdf(file=paste("test","pdf", sep="."))
plotTracks(track)
dev.off()

test.pdf looks like this:

However, what I want is for the boxes to be coloured rather than labelled with the colour word, how to do this?

Note that my real example is much more complex than this.

Did you try the fill argument?

Something like this:

track <- AnnotationTrack(start=c(1,5,7), end=c(2,6,10), strand=c('*','*','*'), stacking="dense", showFeatureId=TRUE, id=c('red','blue', 'red'), fill=c('red','blue', 'red'))