How does sequence alignment help with protein folding prediction? I've been seeing using sequence alignment to help with prediction of protein folding, but am unable to understand why. Google returns other research papers for these search terms, hence I'm hoping the stackoverflow community can help thank you.

Protein sequences from the same family form a sequence alignment. Some positions in the alignment will be conserved, i.e. the same in all sequences. These are often associated with function, e.g. a catalytic residue at position 50.
However there is a more subtle signal in sequence alignments called covariation. Two positions, say residue 30 and residue 60, may appear to change together across the alignment. We infer that they are close in 3D space in the structure:

Initially statistical techniques (DCA, PSICOV etc.) were used to separate direct from indirect (chained) couplings and extract predicted contacts from sequence alignments. Now deep neural networks have improved performance and even allowed the prediction of residue-residue distances and orientations.
These constraints on the structure allow model building, for example using the CNS software. A good set of constraints effectively guarantees the correct fold. The depth and quality of the sequence alignment is the major factor in the accuracy of the predicted constraints.
This approach to protein structure prediction is the current state of the art. Use of covariation from sequence alignments is the only approach that has reliably given accurate models for larger proteins without templates. See this paper and this paper for the best current methods, and tune into the CASP conference this December to see if anyone can improve on them.


How do you query and explore ENCODE data? I am looking for a modular way to query data from ENCODE.

For example, I would like to get CHiP-seq or similar tracks for a specific cell line. What's the proper way to do it?

Finally, is there an API to do it?

I don't know whether there is an API, but ENCODE's website does provide an interactive data matrix where you can filter data based on assay and sample type, place data sets in a "shopping cart", and then proceed to "checkout" to download the files of interest.


Are there any resources for calculating / representing evolutionary relationships within gene clusters? I am working with a gene cluster found in several species. The genes involved are pretty much the same, with the odd insertion, but they can be often found in different orientations. So, I can't just do a multiple alignment of the whole cluster. I've done phylogenies of individual genes and I've looked at cluster architecture (indels, rearrangements) and things make sense (eg. limited HGT), but..

...Is there a simpler way to do this? A way to calculate / represent a phylogeny of the whole cluster? Like...some sort of web, maybe? Most of the software that does take rearrangements into consideration is meant for whole genomes, and is unreliable when it comes to fragments about 20-30kb long (like mine).

Thank you for your help!

Firstly you wanted a webserver to build a tree, ok, I would use PhyML server. The good thing is that with only 22 taxa the calculations are easily doable.

What I'm assuming is you can map your rearrangements to produce a single homologous gene alignment for all taxa. So you appear to infer you can break your 20-30kb contig into individual aligments of homologous genes which has incorporated the rearrangement and enables easy production of an alignment file. This would leave you with around 20-30 gene alignment files. I hope that is correct, because that is how you get around the rearrangement problem and enables the data to be reassembled for the final tree (below).

Doing the tree is simple, because you can exploit the assumptions used in tree building regarding independence between mutations at each alignment position. The theoretical proof, the data is sound, is complicated. In summary,

build a consensus phylogeny (one tree representative of lots of trees) using two different approaches
combine the data into a single alignment, if the above looks good, or excluding the alignments that are affected recombination.

The theoretical issues with what you want to do are complicated because you need to distinguish bifurcating (clonal) homologues from recombination that could have occured in the reassortment events. The following is a  purists approach, you might not want to do this stuff, but it gives you some idea.

Just before I do that ... I just want to check the alignments assumption is correct? and which bacteria?


Is there a JSON-based genomic feature format? There's a variety of formats for storing genomic features (e.g. Genbank, GFF, GTF, Bedfile) but all of the ones that I'm familiar with use either a custom format (Genbank) or a CSV with defined column conventions and space for additional data (GFF, GTF, Bedfile).

Is there a feature format that is JSON based and adheres to a defined schema?

JVCF is a JSON schema for describing genetic variants. BioJSON is a JSON schema for multiple sequence alignments. The output of mygene.info's gene annotation web service is JSON-formatted.

You might build one of your own schemas from any or all of these, as inspiration for structure and content.

Consider that most genomics formats are historically tabular, because they are not hierarchical, because they originated prior to the advent of JavaScript as a common programming language (and subsequent use of JSON to serialize data), and because a lot of stream-oriented UNIX tools read in data and process it in ways that are natural to tabular structure.

Some tabular results can be hierarchical, but use inefficiencies in field values to apply hierarchy. JSON is inefficient for storing tabular data, but it is absolutely wonderful for hierarchical data. This difference is perhaps worth considering.

One could think of storing a gene as a master object containing various haplotypes, exons, introns, features, associated variants, associated TF binding sites, gene ontologies, etc. Ordering or prioritization of such items would be possible with the use of sorted lists within a gene object. The format is as extensible as it is a matter of adding a property, list, or object to the schema.

To gain adoption or "buy-in", build tools which publish it to users, and tools that consume what users get in order to do things, along with conversion and validation utilities.


Uppercase vs lowercase letters in reference genome I am using a reference genome for mm10 mouse downloaded from NCBI, and would like to understand in greater detail the difference between lowercase and uppercase letters, which make up roughly equal parts of the genome. I understand that N is used for 'hard masking' (areas in the genome that could not be assembled) and lowercase letters for 'soft masking' in repeat regions.

What does this soft masking actually mean?
How confident can I be about the sequence in these regions?
What does a lowercase n represent?

What does this soft masking actually mean?

A lot of the sequence in genomes are repetitive. Human genome, for example, has (at least) two-third repetitive elements.[1].

These repetitive elements are soft-masked by converting the upper case letters to lower case. An important use-case of these soft-masked bases will be in homology searches: An atatatatatat will tend to appear both in human and mouse genomes but is likely non-homologous.

  How confident can I be about the sequence in these regions?

As you can be about in non soft-masked based positions. Soft-masking is done after determining portions in the genome that are likely repetitive. There is no uncertainty whether a particular base is 'A' or 'G', just that it is part of a repeat and hence should be represented as an 'a'.

  What does a lowercase n represent?

UCSC uses Tandom Repeat Finder and RepeatMasker for soft-masking potential repeats. NCBI most likely uses TANTAN. 'N's represents no sequence information is available for that base. It being replaced by 'n' is likely an artifact of the repeat-masking software where it soft-masks an 'N' by an 'n' to indicate that portion of the genome is likely a repeat too.

[1] http://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1002384


Master of Pores Installation Mac OS Mojave 10.14 I'm bit of a newbie with ONT analyses and am having some trouble with the Master Of Pores installation instructions (https://biocorecrg.github.io/master_of_pores/install.html). I've managed to install nextflow, cloned the repository and have downloaded albacore and guppy. But when I run the pipeline below, i get an error saying
nextflow: command not found. Any help would be most appreciated.
curl https://get.nextflow.io | bash
mv nextflow ~/master_of_pores/NanoPreProcess/bin
nextflow run nanopreprocess.nf

you should not move nextflow within ~/master_of_pores/NanoPreProcess/bin .
Just add nextflow in your path.
export PATH=$PATH:where_is_your_nf_executable

Feel free to open an issue on github if you need more help.


How to learn/get started bioinformatics? I know this is a broad question and I'm sorry if this is not the right place to post (this is my first post ever), but is there some sort of roadmap for bioinformatics?

I'm an undergraduate student and I have a strong background in statistics, computer science, and molecular biology so I feel like I should be able to get around, but every time I talk or read about bioinformatics I feel lost. It's such a vast area that I don't even know where to begin or what even exists.
Is there some kind of all-encompassing textbook or roadmap for me to follow?

Most textbooks I could find expend way to much time explaining really basic stuff like what is a cell and how to program and I feel they fail to convey the bigger picture of bioinformatics. Things like what exists in the area, what are the current technologies, conventions, and challenges that we face. What should I know in order to call myself a bioinformatician?

The formal answer is O'Reilly "Bioinformatics programming, using Python". Check out the classic O'Reilly texts on "Mastering Perl for bioinformatics", and possibly"Beginning Perl for bioinformatics" (which is likely far too basic), which establish some core principles and have a solid reputation in teaching. However, the Perl/Python thing might get confusing (I dunno).

If you are statistics orientated at a guess, O'Reilly's, Bioinformatics data skills... I don't know, but something that would focus on R.

The real goal is core-technical competance because the biological understanding will be taught by your supervisor. Thereafter you might turn your attention to teaming up with a lab for project training.


What unit do I get on the y-axis of a metagene profile plot? I start with a sorted and indexed bam file ("mapped.bam") representing the mapping of small reads on a reference genome, and a bed file ("genes.bed") containing the coordinates of a set of features of interest (let's say they are genes), for which I want to compute an average profile using programs from deeptools. I would like to understand the steps involved to be sure of what the vertical axis of the final profile represents.

First step: making a bigwig file

I create a bigwig file ("mapped.bw") from the bam file using bamCoverage as follows:

bamCoverage -b mapped.bam -bs 10 -of=bigwig -o mapped.bw

The help of bamCoverage says:

  The coverage is calculated as the number of reads per bin, where bins are short consecutive counting windows of a defined size.

In my case, the bins are 10 bp long. My reads are longer than that.

For a given bin, a given read can:

completely overlap the bin
overlap the bin on n bp, n < 10
not overlap the bin at all

Please correct me if I'm wrong: My guess it that the read is counted as 1 in cases 1. and 2., and 0 otherwise, and I also suppose that a read can be counted for several successive bins if it is long enough.

Second step: averaging over genes and plotting

I compute a "meta profile matrix" ("mapped_on_genes.gz") using computeMatrix scale-regions as follows:

computeMatrix scale-regions \
    -S mapped.bw \
    -R genes.bed \
    --upstream 300 \
    --unscaled5prime 500 \
    --regionBodyLength 2000 \
    --unscaled3prime 500 \
    --downstream 300 \
    -out mapped_on_genes.gz

(There is a -bs parameter which default value is 10 according to the help of the command.)

I use this to plot a profile using plotProfile:

plotProfile -m mapped_on_genes.gz \
    -out mapped_on_genes_meta_profile.pdf

I obtain a profile in with values on the y axis. In what units are these values?

My guess is the following:

For the upstream (300 bp) and internal 5-prime (500 bp), since the bin size was the same in bamCoverage and computeMatrix, each point on the x axis probably represents a 10 bp window, and its y coordinate is the average over the regions present in the bed file of the corresponding bins in the bigwig file, so it is an average number of reads overlapping a 10 bp bin.

Same thing at the 3-prime and downstream side.

For the central 100 bp portion, before averaging over regions some shrinking or spreading of the bins must have been performed, I guess by averaging between neighbouring bins. So the final unit is still a number of reads overlapping a 10 bp bin.

And if I use larger bins, I should end up with proportionally higher values.

Am I correct?

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

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

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

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

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


validating identified sub-populations of cells in scRNA-seq In the analyses of single-cell RNA-seq data there are different unsupervised approaches to identify putative subpopulations (e.g. as available with Suerat or SCDE packages).

Is there a good way of computationally validating the cluster solutions? Different methods may results in slightly different clustering results. How to know which one is the best i.e. representative of biological sub-populations?

A SC3, single-cell consensus clustering, approach could be helpful here. It aims at achieving "high accuracy and robustness by combining multiple clustering solutions through a consensus approach" https://www.nature.com/nmeth/journal/v14/n5/full/nmeth.4236.html


Best way to "fish" long reads that match a query sequence Very simple set-up: We have a PacBio long high-fidelity (HiFi) reads, genome sequencing, and some of those that contain a particular sequence. I want to find out which ones. The length of the query sequence is 7kb, though ideally I would use as query a sub-sequence of only 1.5kb.

The output can the just the read ID or the full read sequence.

We tried DAMAPPER and it gives some decent results but it might not be appropriate to search for the shorter sub-sequence of 1.5kb. It also feels like an overkill for this purpose.

Bear in mind that this is my first foray into long-read sequencing so my mindset is still very much one of short, low error-rate reads.

I have seen the blasr aligner typically used for aligning PacBio reads to a reference.  Is there a reason why blasr would not be a good choice? I guess the very short sequence you are trying to align to might be an issue.

Minimap2 is capable of mapping pacbio reads, and can also do pairwise mapping to see if reads overlap with each other, so this function might help you:
From https://github.com/lh3/minimap2

minimap2 -x ava-pb  reads.fq query_seq.fq > ovlp.paf


PacBio long-reads impact in transcriptome de novo assembly? We are strongly interested in assembly a good transcriptome of reference for a non-model organism and build a local database. We have sequenced the same individual with Illumina (150 millions of pair-end reads) and PacBios IsoSeq v3 (2 SMRT cells, one for shorter transcripts, shorter than 5kb and other for longer transcripts, up to 5kb).
To process long-reads, I have followed the PacBio IsoSeq pipeline proposed in their Github repo (https://github.com/PacificBiosciences/IsoSeq). The final result was removing 70% of the long-reads. Is that normal?
Using this data, I have assembled the transcriptome using only short reads and another combining long- and short-reads. In the end, I have not found any difference... Approx. The same N50, the same number of transcripts assembled, the rate of misassemblies... Does anyone know if PacBio data does not worth for transcriptome de novo assembly?

A few comments:
Never use N50 as a metric especially for transcriptomes. It has some semblance of relevance for genome assembly, but all that is void for a transcriptome with inherently dynamic lengths.
At the end of your IsoSeq pipeline, you should have (ideally) full-length transcripts. Have you considered that you don't miss much in your IsoSeq data? The biggest gain using IsoSeq is that it should have fewer non-existent isoforms assembled than a De-Bruin graph approach.
Let's say you have the following locus in your genome:
5'                                                                 3'
------------=======-----------=======----------=======-------------
Upstream      E1     Intron     E2    Intron     E3      Downstream

There are theoretically 7 isoforms:

E1
E2
E3
E1,E2
E2,E3
E1,E3
E1,E2,E3

A short read assembler has to look at evidence that connects exons to see which isoforms are represented. Let's say you have the following PE reads:
5'                                                                 3'
------------=======-----------=======----------=======-------------
               -->..............<--- Read 1
                                 -->...............<-- Read 2
                 --___________-->..<-__________--- Read 3

Read 1 and 2 both have a mate in either exon, so that's good evidence that they come from a single transcript. You can also have reads like Read 3, which spans two exon-exon junctions, again good evidence. There are some problems, however, if things get more complicated: if the cumulative length of your exons is greater than the fragment length of your short reads, or you have many exons in a gene, you will never get a single read pair that confirms the existence of some transcripts. The same goes for the opposite. If you find a read pair on a single exon, that is not evidence that this single exon is a transcript in and of itself.
Here's where short read assemblers take different approaches. Some make educated guesses about the existence of some transcripts, e.g. by utilizing expression counts in an EM or dynamic programming approach, some just enumerate possibilities.
And that's where IsoSeq comes in. A single, long read from a single transcript helps a lot to add evidence and make sure your transcripts actually represent the transcripts that exist in your organism. So, often, your IsoSeq should be fewer transcripts than your short read assembly.
In essence, what you should do, is take the IsoSeq and curate those as a set of high-confidence transcripts, i.e. those that you have observed from a single read source. That's not to say that some transcripts that exist only in the short read assembly aren't real, just that you have less definitive evidence for them.
If you are mostly interested in the gene level, that doesn't make a huge difference, but once you need to look into transcript level usage/expression etc it helps a lot to be able to eliminate low confidence transcripts.
Use tools better suited for transcriptome assembly evaluation:
Detonate, TransRate and (Busco) and compare your assemblies using those.


FASTQC overrepresented sequences after trimming I have a set of RNA-seq samples from different experiments (Single and Paired End, depending on the experiment). I ran FASTQC in all the samples and found overrepresented adapter sequences:

I removed the adaptors (TruSeq adaptors) using Cutadapt (in addition, I removed low quality and N bases from the 3' end of the reads). After that, I ran again FASTQC and the output is the following (representative example) :

Does anyone know what is happening? Now I have an overrepresented sequence for which no sequence is provided. What does this mean?

As @AaronBerlin mentioned, you didn't remove reads that were completely trimmed. Next time use the --minimum-length option and set it to something reasonable, like 20. Alternatively, use "Trim Galore!", which is a wrapper around cutadapt that has more reasonable defaults.


Determining the quality of a publicly available genome assembly I am looking for homologous genes in various protozoan species using BLAST. The genome sequences of these species are deposited in NCBI's WGS database. NCBI’s webpage includes the global statistics for each WGS project and it indicates the assembly level (scaffold, contig, etc.), the assembly statistics (number of contigs, N50, L50, etc) and (if available) the genome coverage. Can any / some of these parameters be used to evaluate the genome assembly’s "quality", especially if there is more than one genome assembly for the same species? Can these data be used to determine if the genome is very fragmented or not?

I understand that quality is a relative term and that there are several options for assessing genome assemblies (QUAST, BUSCO, etc.) but I am just looking for a "simple" way to compare several genome assemblies (for e.g., this post How can I improve a long-read assembly with a repetitive genome? states that the reference genome is highly fragmented, but there is no explanation regarding how the OP determined this). I believe that the N50 value can be helpful, but I need some context to understand what is "good" or "bad".  I am at a loss of where to start, so any information will be greatly appreciated!

-Leah

It is not so easy to determine the quality of an assembly - but sequence length, number of contigs, N50 and L50 can give you some guidance.

How much does the sequence length deviate from the expected length (or how much variation is among assemblies of the same species). Same for number of contigs/scaffolds. The sequence length, and number of contigs together with N50 and L50 give you an idea about the fragmentation. Very small N50 compared to expected genome size and chromosomes is a likely indication of high fragmentation, same for high L50 compared to expected number of chromosomes.

However, these number can be misleading. High N50 and low L50 can occur witho some very long contigs, so you may want to compare N75/L75 or better the Quast statistics


Classification (supervised learning) of expression data on pathway level I was wondering if there is any way to apply classification algorithms (e.g random forest) on microarray data but not using the genes as predictors/features but the pathway they belong to.

The thing is that the expressions of the genes that belong to the same pathway should somehow, grouped together and transformed into a new value that will describe that pathway for each specific sample.

Let me give an example.

Assuming that we have g1, g2, g3, g4, g5, g6, g7, g8, g9, g10 as genes
and p1, p2, p3 as pathways.

And g1, g2, g3 belong to p1 while g3, g4, g5, g6belong to p2 and g2, g6, g7, g8, g9, g10 belong to p3

One thought is to calculate the mean expression of the above gene groups per sample and then run the classification algorithm on top of these values. But I would like to hear/read your opinion on such an approach and also your suggestions on using a different method to do this.

The final goal would be the extraction of the most significant classification predictors (pathways).

Have a look at the GSVA package. It allows to convert a matrix with genes x Samples to a pathways x Samples using several methods ssgsea, gage, gsva...

Afterwards you can use that matrix as input for differential expression of pathways or classification algorithms or whatever.

However it depends on the input of the "pathways" you give it. Make sure that the pathways are correct. Also consider that the pathway topology might change with conditions.


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

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

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

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

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


Brackets in .Rmd code blocks I am making some R markdown documents for my stats class.  I've noticed I get an error in some functions when I have a code block with {r}  instead of just r.  Like so...
``` r
this code will be fine
```

``` {r}
this code will have an error
```

When should I be using the brackets for codeblocks in .Rmd and when should I leave them out?

Does it work without a space between the ``` and the {r}? It may be that the Rmd processor is getting confused about the spaces, but ignoring the bare r, and assuming the code is R by default.
The examples here don't have spaces:
https://rmarkdown.rstudio.com/lesson-3.html


Human Genome Alignment and Variant Calling Benchmarks I am doing a project that requires bench-marking the performance of computational tools for alignment and variant calling of human genome sequencing data. In particular I would like to establish benchmarks comparing the performance of GATK and ADAM/Avocado.

My main question is which genomic data-sets should I use as the input to the data processing pipeline. In Frank Nothaft's (author of ADAM and Avodaco) thesis he reports bench-marking the processing of NA12878 from the 1,000 genome's project in Avocado as compared to GATK's HalotypeCaller. He also reports benchmarking on 270 samples from the Simons Genome Diversity Data-set.

I am interested in reproducing this portion of his work, but do not know where I can get a hold of these datasets. Where can I get a copy of the NA12878 genome (BAM file) and also were can I get aligned samples from Simons Genome Diversity Data-set?

NA12878

My recommendation: download raw reads from Illumina BaseSpace and do alignment yourself. Aligning ~40X worth of human reads takes overnight with ~16 CPU cores.
You can acquire Platinum Genome data from ENA. There is a BAM available for download, but that was done by Illumina's own aligner against a wrong reference genome. You'd better do alignment by yourself. Another problem with this dataset is it is too good to be representative. Modern data are worse.
GIAB provides BAMs here:

ftp://ftp-trace.ncbi.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/

It is not "typical", either.

SGDP

Instructions on downloading SGDP BAMs are available on this page, in particular this Word document. Warning: you need ~30TB disk space for these BAMs. I also heard it is tricky to download data.
There is apparently a mirror from cancer genomics cloud. I don't know how easy it is to access.


What are the standard ways to visualize protein-protein or gene-gene interactions? I would like to visualize this interaction list.  Is there an online/web based way to do it? Is there a way to analyze the data?

Assuming I have exported the list into lines such as:

2  A2M 348 APOE

since I am not yet allowed to comment, I will have to pots this as an answer.

I guess Cytoscape (http://cytoscape.org/) would be a nice way to analyse such interaction lists/networks (still depends a bit on what you plan to do).

However, for this you will have to preprocess the downloadable excel sheet a bit.

save as tab-separated file (tsv)
run the below python3 script to format into an easily importable format for cytoscape. You can run the script as python3 script.py <path-to-tsv-file> <delimeter for cytoscape> . The new format is written to std out on the console, therefore maybe you need to write the output to a new file. A sample call could look like python3 testbench.py human_interactome_2017-06-12.tsv.csv ";". Make sure to add the delimeter in quotes.
you can then import the new file into cytoscape (import network as table) where the first column is your source, the third column your target and the middle column your interaction type.
Within cytoscape you can then adapt your edge style to vary the line type in a discrete mapping depending on the value of the interaction.

Screenshot from cytoscape:

Code below:

import sys

HIIIidx = 4
HII5idx = 5
venkatesanidx = 6
yuidx = 7
Litidx = 8

print(sys.argv)

if len(sys.argv) <= 1:
    print("Usage: script.py file delimeter")

delim = ';'
inFile = sys.argv[1]

if len(sys.argv) > 2:
    delim = sys.argv[2]

def printInteraction(syma, symb, inter):
    print(syma + delim + inter + delim + symb)

with open( inFile , 'r') as file:

    file.readline()

    for line in file:

        aline = line.strip('\n').split('\t')

        if aline[HIIIidx] == '1':
            printInteraction(aline[1], aline[3], 'HIII')

        if aline[HII5idx] == '1':
            printInteraction(aline[1], aline[3], 'HII5')

        if aline[venkatesanidx] == '1':
            printInteraction(aline[1], aline[3], 'VEN')

        if aline[yuidx] == '1':
            printInteraction(aline[1], aline[3], 'YU')

        if aline[Litidx] == '1':
            printInteraction(aline[1], aline[3], 'LIT')


What is the detection limit of protocol in scRNA-seq In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

  Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA.

I used to think that the detection limit of a protocol was the average maximum number of molecules detected across cells in a scRNA-seq experiment. However, by working on the dataset provided by The Tabula Muris Consortium et al. (pp) I noticed that RNA raw counts often spikes to more than 10 molecules.

What is the actual meaning of detection limit of a protocol in scRNA-seq?

The detection limit refers to the ability to detect a given gene/transcript. In this case, they're saying that a given transcript requires 5-10 copies in order to be likely detected (in other words, the loss rate is 80-90%, which is directly related to a previous question of yours).


Pipeline to analyse small RNAseq and tRNAs? I am a newbie in this world (just finished my master's on bioinformatics. As a practice, a researcher gave me 4 fastq files from a small RNAseq experiment, to see if I am able to reproduce their results but with an updated pipeline and tools (they used bowtie in 2014).
I used fastqc for the quality check and I removed the adapters properly.

My first intention was to use Salmon for pseudo-align them, but after do some reading it seems not to be a good idea for such small RNAs.

Can anyone provide me with some light on this? Which alignment tool use and which parameters, and etc. Ideally a proper protocol/pipeline would be amazing.

Thank you in advance...

PD: I was thinking to use the HISAT2 with the full genome (or the pre-builded index from HISAT2 website).

Kind regards,
and thank you in advance!

I have heard that the Small RNA Workbench from UEA can be an easy solution for a full small RNA pipeline.  It is java based, has a GUI, and is supposed to be easy to install and run.  I am not an expert in small RNA analysis, but I think this pipeline will certainly given you an idea about the steps needed for a small RNA analysis, it has been put together by experts.


Negative scale-free topology I am trying to build gene co-expression networks with WGCNA by combining in-house and publicly available RNA-seq datasets. I am interested in identifying gene networks associated with different environmental conditions. I used SVA to adjust for batch effects and "protected" different environmental conditions in the model as I don't want to lose information about them. However, when I used the "batch-adjusted" data to create networks, all the scale-free topology indices were negative (max = -0.05).

I tried generating the networks with softpower = 16 as it is suggested for sample sizes 20-30. However, only five networks were generated which is far less than I expected.

What might be the reason for getting negative scale-free toplogy indices?

Can combining studies from different environmental conditions cause the networks to be not biologically meaningful?

Whenever you have batch or other nuisance variables driving variation you're going to get odd networks. You tried using SVA, which is a good method to use, but it seems in this case that it's able to optimally handle this dataset. I also strongly suggest that you use pickSoftThreshold() to choose a softpower, since that may produce better results (though since you already have negative values I'm not holding by breath). In general, you're going to need to play with either SVA or filtering to remove the batch-related genes that are causing problems.


Sequence alignment using BWT My Problem:
Skipping some specific background, what I want to do is judging whether some soft-clipping sequences are the same, which may result by the same SV event. Colored bases in Fig.1 is an example of soft-clipping sequences.

My solution:
I suppose this is a sequence alignment problem. Then I use FM-indexing based on BWT to do this. My inexact/approximate matching implementation is referred from BWA.
My workflow:

Extract these soft-clipping sequences
Choose the longest one as the reference
Iterate others to match with the reference, if the number of unmatched is large than 1, it means these sequences are not the same.

Question 1:
Is that sounds reasonable to you?  Any suggestions?

Question 2:
If it is OK, how should I implement the penalty scoring scheme?
Here, because these sequences should come from the same SV event, the main differences(Mismatch, InDel) should arise from sequencing errors, not small genomic variants. So I want to set a higher gap open penalty and even higher gap extend penalty, like
mismatch_penalty = 4
gapopen_penalty = 11
gapextend_penalty = 20
alignment penalty cutoff = 30

The later means that one alignment with two more gaps will be considered as spurious alignment.
How do I optimize these parameters?

You just need to check the flag for suplementary mapping. If there is suplementary mapping then use the coordinates; if they are all the same the seq is the same and you just have one event, if different reads have different suplementary mapping locations then there's multiple events. The problem is more complex than this though so just use an existing tool.


How to interpret PCA output statistically and biologically? How can I interpret the PCA results statistically for biological data?

I have used FactoMineR and factoextra libraries for PCA

Scripts used:

library(FactoMineR)

res.PCA = PCA(df, scale.unit=TRUE, ncp=4, graph=F )
par(mfrow=c(1,2))
plot.PCA(res.PCA, axes=c(1, 2), choix="ind")
plot.PCA(res.PCA, axes=c(1, 2), choix="var")
dimdesc(res.PCA, axes=c(1,2))

library("factoextra")
fviz_pca_var(res.PCA, arrowsize = 1, labelsize = 3, repel = TRUE, col.var = "contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

Here is my PCA output of differnt KO genotypes); how can I explain/interpret the PCA output statistically and if possible biologically: esp, how this colors /contrib explain

All that plots like this are telling you is that there are some genes that contribute more to the variability seen between your various samples than others. In an ideal world these genes will also be differentially expressed between your groups, but since we don't live in an ideal world that might not be the case. In general, the longer the line and the closer it is to the X axis the more it contributes to the overall variation (since it's contributing more the PC1, which is 53.5% of the variation in your data).

Since you've asked in the comments for a brief "results and description" section style summary of such a PCA I'll include something brief below. Note that one should not over-interpret PCA plots. The axes are (mathematically) orthogonal and don't necessarily correspond to anything biologically relevant.

  We performed a PCA on the variance-stabilized counts to check for batch effects and overall clustering of the data. As can be seen, the "3T" and "5T" groups cluster together along the first principal component, while the "0T" and "1T" samples cluster on the opposite side. This suggests that these groups are more similar to each other, as further evidenced by ...

Assuming you had used the vst() function in DESeq2 (I note you've tagged your post with RNAseq, so perhaps you're using DESeq2) then something like that above would be a reasonable description.


chipseeker annotation issue library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)

peaks <- readPeakFile("MACS_folder/pol2_chip/Pol_2pausing/pol2_diff_peak.bed")

print(peaks)

#peakAnno <- annotatePeak("merge/blast.bed", tssRegion=c(-3000, 3000),
#                        TxDb=txdb, annoDb="org.Hs.eg.db")

peakAnno <- annotatePeak(peaks, tssRegion=c(-1000, 1000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")

So this is basically what im doing to annotate the peaks.
Now the issue comes with annotation for example I see TLR4 annotated to this
118082896   118083202   307 *   shp1_R1_peak_679    Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9   117708197   117708979   783 1   7099    ENST00000490685.1   374699  ENSG00000136869 TLR4

Now The peak is 679 as you can see from the IGV plot as well as annotated file. But when i see in IGV browser TLR4 is far away from the peak where i get the annotation.
So now is it an issue with the code im using or due to the inherent nature of the data.
Any help or suggestion would be really appreciated.

I had a similar case just recently.
The tssRegion argument ChIPseeker is used to define what is considered a promoter, which is frequently the question at hand. However, ChIPseeker will not perform any kind of filtering on the annotated peak file.
In the results there is a column annotation that contains information what kind of genomic region the peak is in. If it is within your tssRegion window it will be labeled as a promoter.
In your case the annotation is 1st Exon = Exon (ENST00000413759.1/ENST00000413759.1, exon 1 of 1) 9. This tells you that the peak is not part of the promoter region you defined and located within the first exon of the indicated transcript. The peak is really far away from TLR4 (374699 bp) but I think all other sequences nearby are pseudogenes and this is why TLR4 is selected as the closest gene.
You can get an overview over the annotations using plotAnnoPie(peakAnno)
tl;dr: you need to filter the results of the ChIPseeker annotation to get the genomic regions you want.


GEOquery errors due to list input I want to extract the expression matrix from the GSE object in GEOquery.

library(GEOquery)

gse1 <- getGEO("GSE97168",GSEMatrix=TRUE)
datExpr = exprs (gse1[[1]])
datExpr

GSM2553017 GSM2553018 GSM2553019 GSM2553020 GSM2553021 GSM2553022 GSM2553023 GSM2553024 GSM2553025
     GSM2553026

Table(GSMList(gse1)[[1]])[1:5,]

Table(GSMList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GSMList’ for signature ‘"list"’

Table(GPLList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GPLList’ for signature ‘"list"’

Meta(gse1)$description
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Meta’ for signature ‘"list"’

Columns(gse1)
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Columns’ for signature ‘"list"’

Why it's not working with R version 3.5.1?

Use GSEMatrix=FALSE, the information you want appears to be in the soft file, rather than the matrix. Yes, this is moderately annoying since it means you need two objects to describe this particular dataset.


Nextflow: ordering the execution of processes; outputs of process A are inputs of process B I need to build a pipeline where 2nd process starts after completion of the 1st process and the 3rd after 2nd. As I understand nextflow's processes run parallel. But I need it to run in a order 1st>2nd>3rd processes and so on. Because each process runs on an array and produce multiple outputs which are used by the next. If they run paralleled, the 2nd process may look for input 20 whereas the 1st process produced output 10 and then the process stops, because the command in the 2nd process does not find the input. Is there any way to accomplish this. I tried as follows:
process A{
input:
tuple val(some_name), path(some path) from some_channel
output:
tuple val(some val), path(some path) into A_some_ch
tuple va(some other va), path(some other path) into A_some_other_ch
val true into done_ch
script:
"""
some command out some files
some other command out some other files
"""
}
process B {
input:
val flag from done_ch
val (some val), path(some path) from A_some_ch
val (some other val), path(some other path) from A_some_other_ch
output:
some val and path
script:
"""
some commands out some files
"""
}

I got the idea using of flag from here: https://github.com/nextflow-io/patterns/blob/master/docs/mock-dependency.adoc
But if I use the flag, I got:
Operation 'set' does not allow more than one target name

Does it clarify the issue? Any Help?
Best, Zillur

Nextflow uses channels (asynchronous FIFO queues) to communicate between processes. The interaction between these processes, and ultimately the pipeline execution flow itself, is implicitly defined by these input and output declarations.
The issue with your code, I think, is the use of multiple input channels. Consider the following example that uses only a single queue channel. To avoid declaring the inputs again in process A's output declaration, we can output a key to join on as the first element in the output tuple. In the case below, this is simply the 'some_name' variable. We can then join up the two channels using the join operator as follows:
Channel
    .fromFilePairs( './data/*.txt', size:1 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .into { input_ch_A; input_ch_B }

process A {

    input:
    tuple val(some_name), path(some_path) from input_ch_A

    output:
    tuple val(some_name), val("procA_${some_name}"), path("procA_${some_name}") \
        into procA_results

    """
    touch "procA_${some_name}"
    """
}

process B {

    input:
    tuple val(some_name), path(some_path), val(some_other_val), path(some_other_path) \
        from input_ch_B.join(procA_results)

    """
    true
    """
}
```


Error: Insufficient values in manual scale. 24 needed but only 1 provided? folks.

Can you give me idea about an issue that I encountered?
I ran single-cell RNA sequencing data using Suerat in R, and when I tried to draw violin plot, there is an error as below.

  VlnPlot(tumors, "Cd79a", "Foxp3", "CD14")

Error: Insufficient values in manual scale. 24 needed but only 1 provided.

Run rlang::last_error() to see where the error occurred.

However, I can draw the violin plot for Cd79a. The problem happened to only several genes. I haven't solved this problem yet. Do you know what should we do?

Thank you.

You must supply genes (or features) as a vector, not individually:

VlnPlot(tumors, c("Cd79a", "Foxp3", "CD14"))

By the way it is probably Cd14, not CD14.


Why would someone use a CRAM instead of a BAM? I had this question from a graduate student yesterday, and I was stuck.

What should I say? Why use a CRAM instead of a BAM?

When is it a good idea to use a CRAM instead of a BAM?
When is it a bad idea?

Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings.
If you need to guarantee that any random obscure downstream program will be able to handle it.

Uptake of CRAM has been pretty slow. Java programs using htsjdk (e.g., picard, IGV and GATK) have only relatively recently added support for CRAM. If you need to use an old version of those for some very odd reason then CRAM may not be supported.

There are a lot of programs written in python that use pysam to open BAM files and these should, theoretically, support CRAM. The issue is that some of the functions may fail and one can't assume that authors will have always written the code needed to handle this. I'll use deepTools as an example, since I'm one of its developers. One of the things about CRAM files is that they (by default) are made such that you require a reference genome in order to construct the sequence field in each alignment. This works fine if you're using a standard genome (htslib, via pysam, can fetch many standard genomes from the web automatically), but if you're not, then you need to specify a fasta file to use for decompression. Every tool, then, needs to add an option for this. With pysam 0.14 and htslib 1.7 this can be circumvented by not decompressing the sequence, but behavior has to be explicitly requested.

Another issue is that many tools will use features from the file index, such as the .mapped accessor, to get the number of mapped reads in a file. CRAM files contain very very little information, so this then fails. Consequently, tool authors need to check for CRAM files and both derive and propagate this information through their functions if it's needed. This can be a time-consuming task (e.g., it took me a couple days to get this implemented in deepTools). Relatedly, samtools idxstats is useless on CRAM files, since there are no statistics stored in the index.

That having been said, it's likely that CRAMs slowly gaining acceptance will eventually make it the standard. It's already a convenient archival format, it's just a matter of time before users can assume that most analysis programs are written to handle it.


Pausing index calculation I have Pol_II data for Ctrl vs Knockout condition with two replicates each. I want to calculate the pausing index for the Ctrl vs KO condition.

Here I have given tag directory as input to Homer.Now what I see is as such,for each biological replicate it calculates the pausing ratio.How can I calculate as a whole like 2 Ctrl vs 2 KO replicate?

I'm giving a sample of the output.

Any help or suggestion would be really appreciated.

dput(head(CTRL_PTPEM_Pol_II_pausing))
structure(list(`Transcript/RepeatID (cmd=analyzeRepeats.pl rna hg38 -strand both -count pausing -condenseGenes -d Control_Pol_II_new_set/ Control_Pol_II_set2/ PTPN6_KO_Pol_II_PTPEN/ PTPN6_KO_Pol_II_set2/ -rpkm)` = c("NM_173803",
"NM_001301039", "NR_132775", "NM_022114", "NM_001011666", "NR_160728"
), chr = c("chr16", "chr1", "chr3", "chr1", "chr7", "chr19"),
    start = c(15395704, 150272234, 13618331, 3069153, 28686083,
    22527485), end = c(15400754, 150277284, 13623381, 3074203,
    28691133, 22532535), strand = c("+", "+", "+", "+", "+",
    "-"), Length = c(250, 250, 250, 250, 250, 250), Copies = c(2,
    1, 1, 1, 1, 1), `Annotation/Divergence` = c("MPV17L|M-LPH|MLPH1|MLPH2|MPV17L1|-|16p13.11|protein-coding",
    "C1orf54|-|-|1q21.2|protein-coding", "SNORA93|-|-|3p25.1|snoRNA",
    "PRDM16|CMD1LL|KMT8F|LVNC8|MEL1|PFM13|-|1p36.32|protein-coding",
    "CREB5|CRE-BPA|CREB-5|CREBPA|-|7p15.1|protein-coding", "LOC105376917|-|-|19p12|ncRNA"
    ), `Control_Pol_II_new_set/ Pausing Ratio` = c(0.982, 1.832,
    1.114, 1.247, 3.476, 0.632), `Control_Pol_II_new_set/ Promoter Reads (-50,200)` = c(0.125,
    0.748, 0.249, 0.873, 0.499, 0), `Control_Pol_II_new_set/ GeneBody Reads (200,5000)` = c(0.13,
    0.338, 0.208, 0.669, 0.032, 0.091), `Control_Pol_II_set2/ Pausing Ratio` = c(0.686,
    1.131, 1.054, 0.816, 2.205, 1.005), `Control_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.1,
    0.7, 0.1, 0.3, 0.3, 0.1), `Control_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.203,
    0.604, 0.089, 0.396, 0.068, 0.099), `PTPN6_KO_Pol_II_PTPEN/ Pausing Ratio` = c(1.08,
    1.362, 1.108, 1.224, 1.943, 1.2), `PTPN6_KO_Pol_II_PTPEN/ Promoter Reads (-50,200)` = c(0.151,
    0.602, 0.151, 0.301, 0.452, 0.151), `PTPN6_KO_Pol_II_PTPEN/ GeneBody Reads (200,5000)` = c(0.125,
    0.392, 0.118, 0.212, 0.141, 0.094), `PTPN6_KO_Pol_II_set2/ Pausing Ratio` = c(1.652,
    0.909, 1.551, 1.395, 1.774, 0.584), `PTPN6_KO_Pol_II_set2/ Promoter Reads (-50,200)` = c(0.429,
    0.536, 0.429, 0.536, 0.321, 0.107), `PTPN6_KO_Pol_II_set2/ GeneBody Reads (200,5000)` = c(0.206,
    0.603, 0.229, 0.346, 0.123, 0.279)), class = c("tbl_df",
"tbl", "data.frame"), row.names = c(NA, -6L))

but I would say the pausing index is something that is invented, commonly used, but with no clear indication of what is the expected distribution or how it can be subjected to test like you've mentioned.

I assume you would like to get your work published, so going by the papers with IF > 5 such as this, and this, and this, the clear strategy is use a boxplot to show the difference between conditions, and apply a statistical test such as mann whitney u or t.test to say, there is a difference.

Let's try this, sorry not very good in R,using your dataframe which i named x:

library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)

df = x %>% select(contains("pausing") | contains("Transcript")) %>%
rename_at(vars(1),substr,1,10) %>%
pivot_longer(-Transcript,names_to = "condition") %>%
mutate(condition=case_when(
grepl("Control",condition) ~ "Control",
grepl("PTPN6",condition) ~ "PTPN6_KO"
)) %>%
group_by(Transcript,condition) %>% summarize_all(mean)

In the data frame above, I have averaged the pausing index between replicates (not very good but then again, most papers have only 1 replicate), and group them by condition.

# A tibble: 12 x 3
# Groups:   Transcript [6]
   Transcript   condition value
   <chr>        <chr>     <dbl>
 1 NM_001011666 Control   2.84
 2 NM_001011666 PTPN6_KO  1.86
 3 NM_001301039 Control   1.48
 4 NM_001301039 PTPN6_KO  1.14
 5 NM_022114    Control   1.03
 6 NM_022114    PTPN6_KO  1.31
 7 NM_173803    Control   0.834
 8 NM_173803    PTPN6_KO  1.37
 9 NR_132775    Control   1.08
10 NR_132775    PTPN6_KO  1.33
11 NR_160728    Control   0.818
12 NR_160728    PTPN6_KO  0.892

Then we plot:

p <- ggboxplot(df, x = "condition", y = "value",
                color = "condition", palette =c("#c72e29","#016392","#be9c2e","#098154"),
                add = "jitter")
p+stat_compare_means(comparisons=list( c("Control", "PTPN6_KO")))


Adding legend to Botton and top annotation in complexheatmap I have a heat map but top and bottom annotation don't have legend
Top_ha = HeatmapAnnotation(Group = anno_simple(top$Group,
                                               col =structure(hue_pal()(length(unique(top$Group))),
                                                              names = unique(top$Group)),
                                           border = TRUE,
                                           height = unit(5, "mm")),
                       "Total Mutation Burden"=anno_barplot(clin$TMB, height = unit(15, "mm"), width=unit(30, "mm"),gp = gpar(fill = structure(hue_pal()(length(unique(clin$TMB))),
                                                                                names = unique(clin$TMB)))))

Bottom_ha = HeatmapAnnotation(Response = anno_simple(clin$Response,
                                                 border = T,
                                                 col = structure(c("purple4", "grey"), names = c("Responders", "NonResponders")),
                                                 height = unit(5, "mm"),show_),
                          "TP53 Status" = anno_simple(clin$TP53_status,
                                                          border = TRUE,
                                                          col = structure(c("firebrick", "grey"), names = c("Mutant","Wild type")),
                                                          height = unit(5, "mm")))

Right_Anno <- HeatmapAnnotation(Signature = anno_text(right$Signature), which = "row")

Heatmap(mat,
        bottom_annotation = Bottom_ha,
        top_annotation = Top_ha,
        name = "Exposure",
        right_annotation = Right_Anno,
        show_row_dend = FALSE,
        show_parent_dend_line = FALSE,
        column_split = c(top$Group),
        column_title = c("Mutagenic", "DDR Impaired", "C>A/T Dominant"),
        col = colorRamp2(c( 0, 0.5), c("white", "red")),
        show_column_names = FALSE)

As you are seeing Response, group and total mutation burden don't have any legend and annotation
Please can you help me to put legend for annotation?
I have tried annotation_lable but says Error: either provide levelsor provide namedcolors.
Also show_annotation_name = TRUE,show_legend = TRUE did not work

Answer from @ram-rs, converted from comment:
For the Group, try using it without the anno_simple. Just Top_ha = HeatmapAnnotation(Group=top$Group, "Total Mutation Burden" = ...), nothing else for Group. See if that generates a legend at all. We can customize later.hing else for Group. See if that generates a legend at all. We can customize later.


Can exons be located outside of the coding sequence? I have a gff file like this (I edited the name):

scaffold_x  source  exon    2987526 2987805 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2987526 2987805 .   -   1   name "foobar";proteinId 68892;exonNumber 5
scaffold_x  source  exon    2988610 2988745 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2988610 2988745 .   -   2   name "foobar";proteinId 68892;exonNumber 4
scaffold_x  source  exon    2989157 2989645 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2989157 2989645 .   -   2   name "foobar";proteinId 68892;exonNumber 3
scaffold_x  source  exon    2990903 2991541 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2990903 2991541 .   -   2   name "foobar";proteinId 68892;exonNumber 2
scaffold_x  source  stop_codon  2990903 2990905 .   -   0   name "foobar"
scaffold_x  source  exon    2991587 2992430 .   -   .   name "foobar";transcriptId 68892
scaffold_x  source  CDS     2991587 2992430 .   -   0   name "foobar";proteinId 68892;exonNumber 1
scaffold_x  source  start_codon 2992428 2992430 .   -   0   name "foobar"

So the start_codon is 2992428    2992430 and the stop_codon is 2990903   2990905; I would therefore expect all the exons within this range, however, e.g. exon 5 is found in 2987526   2987805 which is clearly outside this range.

Is there something wrong with the gff file itself or is there a plausible explanation for that?

EDIT:

In the corresponding protein fasta I do find an amino acid sequence associated with this prtein ID and name:

>68892|foobar
MHTGD...

Yes, of course. Exons are not limited to the protein coding regions. Many UTRs are in exons. In fact, you even have various cases of UTRs being multiple exons, and being spliced.

What is strange in your file is not so much that you have exons beyond the stop codon, but that you also have them marked as CDS (coding sequence). That isn't possible, no. While there can indeed be exons in the UTRs, those are not coding and shouldn't be called CDS. That's probably just a minor nomenclature issue though and what you are looking at is a spliced UTR.

So, the exons (or the parts of the exons) that fall between the start and stop codons will be translated, while those that fall outside that regions will not be (UTRs).

For example, this is the exonic structure shown for ENST00000617185, one of the transcripts of the human P53 gene:

The boxes are exons and the lines are introns. The colored boxes are protein coding exons, while the empty ones are UTR exons. Note how they are also spliced and that this particular transcript has 3 non-coding exons on the 5' end and one on the 3' end.


Information about control data I am putting this question because I did not find any useful information from internet because of limited access. My question is related to control (or normal) data that we use for somatic mutation detection (with the help of mutect or SomaticSnipper etc) with tumor data.

I don't have control data and my aim is to find somatic mutations in Multiple myeloma WES caner data.

Here are the questions:

Can we use Hg19 as a control data? If not, what are the technical reasons of not using hg19 as a control data because it hg19 is consensus call.

Is there any way to infer (or predict) somatic mutations (either statistically or by something else) out of all possible mutations? For this reason, I want to know what information is given in control data.

Without control data from your subjects, I don't think there's really no way to distinguish somatic mutations from germ-line mutations. The best you can do is to screen out common variants, which are germ-line mutations that are shared by large numbers of individuals using the population frequencies from something like the  Exome Aggregation Consortium:

annotate_variation.pl -downdb -webfrom annovar -build hg19 exac03 humandb/
annotate_variation.pl -filter -build hg19 -dbtype exac03 example/ex1.avinput humandb/


Download program for ΔG for TM insertion I am trying to calculate the ΔG of TM insertion for many protein sequences.

Has anyone come across the source code or a programmatic version for the von Heijne ΔG for TM insertion? The web server is at http://dgpred.cbr.su.se/index.php?p=instructions

I can't find a link in the papers or on the site.

The application is stored in github (ElofssonLab/dgpred) with a link to the website.

To find it I searched for: dgpred download github


How can I take cell number into account to find total RNA expression? I'm hoping to quantitatively show differences in total RNA expression for a gene in a cluster of interest between different experimental groups. My exported average RNA values for each experimental group in my macrophage cluster are shown below:
00 d Control: 10.43
00 d PLX5622: 14.48
07 d Control: 16.08
07 d PLX5622: 17.12
28 d Control: 15.90
28 d PLX6722: 19.40
However, the violin plot for Cd68 in the macrophage cluster looks like this

VlnPlot(combine.combined, features = ("Cd68"), pt.size = 1, idents = "Macrophages", group.by = "orig.ident", split.by = NULL, assay = "RNA"):

The violin plot suggests that the Control groups actually have more total CD68 RNA in them than in the PLX5622 groups, which isn't what the average RNA values show. Any code you can share showing me how to get total RNA per group would be an amazing help!
Thanks so much!

(I currently cant comment as I used up all my rep for a bounty.)
Could you explain how you calculate the average? In you plot the averages are between an expression level of 2 and 3, while you report values over 10. Given the data you show, there is no significant difference and I would guess p-values are above 0.8!
And you mention that the violin plots suggest the control groups have higher expression.. but how do you see that? For me they are identical with minimal trends that the PLX5622 groups are higher.
And to address your question regarding showing the total: No, you would want to show the mean or median, as the total is pointless, given deviating sample sizes.
EDIT: Overall, I have a feeling you uploaded the wrong picture!
EDIT2: I think I got your point now: The difference in mean is caused by an increased number of outliers in your controls at the value 0. Despite those, the means are identical and this is what can be seen through the violin plots. You should rather try to diminish the effect of these outliers.


Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3 I downloaded the annotation of the C. elegans genome in GFF3 format from Ensembl.
I typed the following command, hoping to get the header of the file (lines starting with #).

grep '^#' Caenorhabditis_elegans.WBcel235.95.gff3

This returned 46798 lines, many many more than the dozen or so header lines. After searching the file for the # character, I noticed many lines that contain only three hashtag (#) symbols.

Why are these lines there?

The three # are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the #.

To grep just the header, which has 1 or 2 # from the line start on, you can use extended regular expression:

$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

Which means: "grep all lines that started with 1 or 2 # and are not followed by #"

If you like to create a new gff3 to get rid of this separator lines you can do this:

$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3

The -v prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")


In GFF3, annotating more than one protein-coding gene (i.e. polycistronic) contained in a eukaryotic mRNA Below is an example of a simple GFF3 file:

1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS 3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS 3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

My question is: if we found another coding sequence (encoding a different protein) range from 3752 to 3904, how should the GFF3 be formatted? It seems to me that the GFF3 file can only allow one protein-coding gene per mRNA. If not, could anyone show me one example? Thank you!

[Edit]:
Sorry, using "CDS" is misleading. I changed the title and contents to "contain more than one protein-coding genes" instead of "more than one CDS". I also change the starting position of the second protein to 3752 (was 3750) to reflect the fact. My sincere apologies.

The GFF3 file format specification doesn't care where annotations are described in the file, although it can make some things easier if the annotations are ordered by the start position:

1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3750    3904    .   -   0   ID=ATRV01011:CDS:1;Parent=ATRV01011.0
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

There are many examples in well-studied genomes of overlapping genes, so the GFF3 format needs to be flexible enough to handle these. Different genes can be encoded across the same region, but on different strands. But different genes can also be encoded on the same strand, sometimes sharing an exon.


Non-exome variants called from whole-exome sequencing data I'm working on rather standard whole-exome sequencing data and treat it the same as whole-genome sequencing data (aligning it to the full GRCh38 reference assembly and calling variants with no exome-related parameters).

I recently noticed that my variant calling results contain numerous variants mapped to a non-exomic region, incidentally, some of them were significantly associated with the phenotype I study.

My questions are as follows:

Are the non-exome variants in my exome data invalid? or could it be that the exome kit captured more than it was designed to capture?
Should I treat raw exome sequencing data any different than genome sequencing data, and if so, what parts of my pipeline should be modified to accommodate exome data?

These are called off-targets due to unspecific binding of DNA to the capturing beads during library preparation. This is normal and expected since capturing is not a perfect process. It is correct to align against the whole genome but you should simply limit variant calling to exons (or basically the on-target regions) and ignore everything else.


Looking for a good tool for repeat masking in long sequences The repeat masker site has a size limit of 90kb.
Where can I mask repeats in a 500kb sequence?

Ideally, I need a thorough repeat masking of known and de-novo found repeats of various complexity. I am working on Windows but can access Unix too.

RepeatMasker can downloaded an run locally. There's then no size limit, though do note that it's quite slow for mammalian-sized sequences (you'll want to split it into chunks and run in parallel in such cases). For a small 500kb sequence that won't pose an issue.


Best way to detect long insertions in bisulfite sequencing data? I am interested in identifying indels in whole genome bisulfite sequencing data (76bp paired end). Currently, I do this by setting the -rfg and -rdg affine gap penalty scores for bowtie2 to more permissive values than the default 5+3N and mapping using the bisulfite sequencing alignment wrapper, bismark.

My question is, what values of -rfg and -rdg will allow me to identify longest possible indels without sacrificing alignment quality? Is it better to set the affine penalty to zero with a high penalty for initially opening an indel (ex. 8+0N)? Or is it better to keep the initial penalty low and having a nonzero penalty for extension (ex. 2+1N)?

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


Sort vcf by contig and position within contig Due to tabix constraints, I need to sort a vcf so that contigs and then positions within contigs occur in numerical order in the vcf.  I don't know if the following will sort positions within contigs, but I tried this:

cat in.vcf | vcf-sort -c > out.vcf

but was informed:

Old version of sort command installed, please run without the -c option.

I've read that can cause problems here. Further, my contigs are't numeric, but are something like this: 'Erumb1_s00000030'.

I see two possible solutions.  (1) I am on Mac OS Mojave, and if there were a way to install the a new unix-style sort, I could use it. I tried homebrew and didn't find anything. (2) there are likely some bash commands that would do it.  If anyone has an idea in either of these directions, I would appreciate it.

You could simply use bcftools sort for it:

$ bcftools sort input.vcf > output.vcf

If you really want to use bash only, you can do this:

$ grep "^#" input.vcf > output.vcf
$ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf

The first command will write the header information to the new vcf file.
The second will sort by contig name and position and append the result to the new vcf file.

fin swimmer


co-occuring genes in my list I have a files which has columns representing genes,sample numbers and sample names (55). The genes are displayed by occurrence in the samples (highest to lowest-55 and 2).I wish to calculate the genes which may co-occur in samples.
I tried doing this with awk 'n=x[$1,$2,$3,$4,$5,$6]{print n"\n"$0;} {x[$1,$2,$3,$4,$5,$6]=$0;}' adding from 1 to 55 columns, but it doesn't work.
An example here:
Nfib    55  A1.S1   A2.S2   A3.S3   A4.S4   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.4   IZ_BH_4.5   IZ_BH_4.6   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Lyn 41  A1.S1   A4.S4   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.9   IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Cdk5rap2    52  A1.S1   A2.S2   A3.S3   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.6   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

The first name is of the gene, followed by number of samples and sample names.
for these 3 genes these are the commonly occuring samples:
A1.S1, A7.S7, A8.S8, A9.S9*, A9.S9, A11.S11, IZ_BH_2.1, IZ_BH_2.2, IZ_BH_2.3,IZ_BH_2.4, IZ_BH_2.5,IZ_BH_2.6,IZ_BH_2.7, IZ_BH_2.9, IZ_BH_2.11,IZ_BH_2.12, IZ_BH_3.2,IZ_BH_3.3,IZ_BH_3.4, IZ_BH_3.5, IZ_BH_3.6,IZ_BH_3.8,
IZ_BH_3.9,IZ_BH_3.10,IZ_BH_3.11,IZ_BH_4.1,IZ_BH_4.2,IZ_BH_4.4,IZ_BH_4.8,IZ_BH_4.11,IZ_BH_5.1,IZ_BH_5.3,IZ_BH_5.4,IZ_BH_5.5,IZ_BH_5.6,IZ_BH_5.7,IZ_BH_5.8

Kindly help.

I think what makes this especially hard is the way your data is formatted.
The problem would be way easier if you had a matrix of presence / absence instead. Each row would be a gene, each column a sample, and the values would be 0 (absent) or 1 (present).
Awk becomes a bit inconvenient for this task, so I'll use python instead. To get such a matrix from your file, you could do the following:
import numpy as np
import pandas as pd

# Get total number of genes / samples
all_samples = set()
all_genes = set()
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    all_genes.add(line[0])
    all_samples.update(line[2:])
all_samples = list(all_samples)
all_genes = list(all_genes)

# Initialize presence matrix with correct dimensions
pres_mat = np.zeros((len(all_genes), len(all_samples)), dtype=int)

# Map each sample/gene to a row/col
row_map = {gene: i for i, gene in enumerate(all_genes)}
col_map = {sample: i for i, sample in enumerate(all_samples)}

# Fill matrix incrementally
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    curr_gene = line[0]
    row = row_map[curr_gene]
    for sample in line[2:]:
      col = col_map[sample]
      pres_mat[row, col] = 1

# Add names and write to file
df = pd.DataFrame(pres_mat, columns=all_samples)
df.insert(0, 'gene', all_genes)
df.to_csv('gene_presence.tsv, sep='\t')

You should get something that looks like:
       gene  IZ_BH_5.5  A8.S8  IZ_BH_2.1  A9.S9*  IZ_BH_5.3  A9.S9  ...  IZ_BH_3.11  IZ_BH_4.7  IZ_BH_4.9  IZ_BH_2.8  IZ_BH_3.5  IZ_BH_2.6  IZ_BH_4.3
0       Lyn          1      1          1       1          1      1  ...           1          1          0          0          1          1          1
1      Nfib          1      1          1       1          1      1  ...           1          1          1          1          1          1          0
2  Cdk5rap2          1      1          1       1          1      1  ...           1          0          1          1          1          1          1

With this format, it is now easy to compute the most frequent samples. E.g. df.sum(axis=0) to get the number of occurence of each sample for these 3 genes. It would also allow you to do different things, such as computing the pairwise similarity between samples.


tabix indexing a GTEx variant file- ':' delimited I have a gtex variant file- the head of which looks as follows:
phenotype_id                                    variant_id
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_13550_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14671_G_C_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_14677_G_A_b38  ...
chr1:15947:16607:clu_36198:ENSG00000227232.5    chr1_16841_G_T_b38  ...

Effectively I would like to be able to lookup the variants inside a particular gene (the file is sorted) and put them in a temporary file:
E.g. If the variants in gene "ENSG00000148481- MINDY3" are on line number 87528225 to 87536766- effectively what I would like is the equivalent to zcat file.gz | sed -n '87528225,87536766p' > MINDY3.txt. However zgrep ENSG00000148481 file.gz is just as fast as the above...
Hence I thought tabix would be the right tool for this-
I would like to tabix index it to make lookups faster. It is gzip compressed and I will firstly have to do:
zcat gtex.txt.gz | bgzip > gtex.txt.bgz

However I am not quite sure how to proceed from there given that the data is not tab-delimited.
As a trial I tried the first 1000 lines:
zcat gtex.txt.gz | head -n 1000 | bgzip > gtex_1000.txt.bgz

./tabix -p bed gtex_1000.gz #index as a bed file
[get_intv] the following line cannot be parsed and skipped: chr1:15947:16607:clu_36198:ENSG00000227232.5        chr1_13550_G_A_b38     ......

./tabix -p vcf gtex_1000.gz #index as a vcf file
Indexing as a bed file results in a warning while indexing as a vcf gives no warning yet either way when I try to retrieve a sequence:

./tabix test.gz chr1:15000:17000

It returns nothing.
I am starting to think that I will just have to write a script that splits on the ':' and writes the data to a new file.... and then index that file- which will take a huge amount of time
... Does anyone know of a trick to index the files with unconventional delimiting?

Tabix files need a name and a genomic coordinate at least. Tabix can be (ab)used by creating a fake genomic coordinate and then it just indexes names but not sure I would recommend this usage. What is your actual expected application of this?


Is it possible to calculate p-values for fold changes of single replicate RNA-seq samples? I'm thinking that this isn't actually possible, but I'd like to check before I write it off.  I have RNA-seq data from cells under three different conditions, there are no replicates for any of the conditions.  The data has been processed with RSEM, and log2 fold changes have been calculated for each control-test pairing using the normalized expected read counts using EBseq.

If possible, I'd like to also calculate the p-value for each of these fold-changes, however, because there are no replicates I don't think that this is possible.  Is there any way to do this?

DESeq2 is able to adjust p-values for single-replicate samples by estimating shared dispersion across all samples. DESeq2 will give a warning, but try its best to carry out the analysis. More information on doing that can be found in the DESeq2 vignette and the DESeq2 manual page.

The results will not be anywhere near the sensitivity that can be achieved when replicate information is available, but it might be good enough if only very obvious large-scale changes are of importance.


Pathway plot in GeneTrail data base So I'm doing a pathway analysis using GeneTrail database. One thing I am not able to get is this similarity measure being calculated and plotted:

I want to know how this is being made, what is the output being plotted because each gene is may be involved in multiple pathways.

I want to be able to make it in R again. Any suggestion or help would be appreciated. How to plot the enrichment result in this form ?

This is a a symmetric matrix with the pathways in both axes. So for each pair a similarity index is calculated (usually it takes into account the number of genes in each pathway and the total number of genes in both pathways). There are several similarity indexes available. A typical one is Jaccard or Sørensen-Dice index, usually 0 means no similarity and 1 equal pathways.

You can do this kind of analysis in R using BioCor (Disclosure: I am the author of this package) which uses the Dice similarity. You will need the name of the genes of each pathway, after calculating the similarity of pathways using mpathSim, you can then plot it with a heatmap (using the heatmap function).


How to select a power for a scale-free topology network In a weighted gene co-expression network analysis (using WGCNA), the soft-threshold power is recommended as a noise filtering. It consists on raising the correlation to a certain number. To decide this power the scale-free topology is estimated for some powers. The function to estimate this prints:

   Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1   0.9300  3.110          0.996  2960.0    3060.0   3970
2      2   0.7510  1.010          0.964  1750.0    1780.0   2900
3      3   0.1730  0.258          0.806  1170.0    1150.0   2280
4      4   0.0942 -0.183          0.713   833.0     782.0   1870
5      5   0.3800 -0.463          0.777   623.0     559.0   1580
6      6   0.5350 -0.656          0.834   481.0     412.0   1360
7      7   0.6270 -0.797          0.872   381.0     312.0   1190
8      8   0.6870 -0.910          0.900   307.0     241.0   1050
9      9   0.7270 -1.000          0.918   252.0     189.0    936
10    10   0.7490 -1.080          0.928   210.0     150.0    841
11    12   0.7850 -1.190          0.948   150.0      98.0    693
12    14   0.8090 -1.280          0.958   111.0      65.9    582
13    16   0.8290 -1.360          0.968    84.0      45.6    497
14    18   0.8410 -1.410          0.973    65.2      32.1    429
15    20   0.8490 -1.450          0.977    51.6      23.0    375

The recommendations of the FAQ indicate that a SFT.R.sq value should be above

  0.8 for reasonable powers (less than 15 for unsigned or signed hybrid networks, and less than 30 for signed networks)

and the mean connectivity below the hundreds.

Others have used the power just as noise filtering without caring much about the scale-free topology fit. I would pick the first power even if the mean connectivity is in the order of thousands because the scale-free topology fit is pretty high, however the slope is puzzling me.

How should the soft-threshold power be selected?

Based on an example question

I would be concerned about the input data that produced these numbers. This dataset has a really odd shape with a dip at the start:

data.df <- read.table("soft_threshold.txt", header=TRUE);
png("scale-free-194.png");
plot(data.df$SFT.R.sq);
dummy <- dev.off();

Here's what I expect this to look like (from the supplementary information in this paper):

It's easier to select the power/threshold when the plot looks like the second situation. With the data as presented, it doesn't look like any value would be suitable. but if I were forced to pick, then I'd choose 12 because it is the first number above the magic threshold (ignoring the weirdness at the start). However, I'll admit that I don't understand the concepts behind the power, and am just following the recommended instructions, so can't explain why it is what it is.


Discordance in gene signature behavior between bulk and single-cell RNASeq The objective of the following analysis is to identify an activation signature of a specific phenotype on bulk RNASeq and to apply it to single-cell RNA-Seq, in order to identify the population of activated cells in scRNA-Seq samples.

In the first step I used DESeq2 to identify genes that increase during different timepoints (I have 4 timepoints):

Bulk samples are normalized using VST (Variance Stabilizing Transformation - DESeq2)
For each timepoint, differentially expressed genes are computed with respect to the previous timepoint (i.e. 2 vs 1, 3 vs 2, 4 vs 3)
Differentially expressed genes that increase with timepoint (log-FC > 1, padj < .05) are selected (n=24)

As a control, I plotted the activation score (computed with scanpy) based on the 24-gene signature applied to the bulk samples from which it has been extracted. As you can see, the score always increases from the first timepoint (blue) to the last (red).

However, when I use the same scoring procedure on a single-cell RNASeq sample containing the same populations, I think something may be wrong.

In the following the violin plot of the extracted activation score on an unsupervised clustering on single-cell data. Excluding cluster 4, which we identified as contamination, while cluster 1 and cluster 0 are supposed to have the correct degree of activation, clusters 2 and 3 have a lower degree of activation with respect to cluster 0 (which we identified as the first timepoint - lowest activation on bulk).

What may be the causes of this behavior? Right now the possible explanation I came out with is that:

Clusters 2 and 3 have global shifts on gene expression that make TPM normalization over-normalize activation genes

Are there other possible explanations for this behavior, and how to check for them?

One explanation could be that your mapping of clusters to timepoints is not accurate. There are other methods you could look at for doing this, for example scMap, scPred, or Seurat v3 (disclosure: I am one of the Seurat developers).

Another possibility is that there is some problem in the calculation of the activation signature. You could instead look at the expression of the 24 genes themselves in each cluster (for example by drawing a heatmap).

The difference between bulk and single-cell here could also be biological. This could happen if the increase in expression of those genes in bulk was due to a change in cell type frequency over time rather than increased expression.


Should PCA be standardized for gene expression? This is a theory/good practice question more than a technical one. If samples are being plotted on a PCA projection of gene expression data, I'm wondering whether it is standard (and if so, why) to center and scale the PCs.

The reason I ask is that in this case the variables (genes) are being measured in the same scale, so I don't know if it would be necessary to center/scale data.

When the gene expression is scaled and centred you reduce the difference between genes.

Imagine you have one gene A that is highly expressed usually and has a standard deviation of 500 units compared to a gene B that is not much expressed and only have a standard deviation of 5.

In the scaled and centred genes both contribute the same because A usually is expressed 10000 and B is usually expressed 100 units. So, for both the standard deviation is 5% of their expression. Meaning that a variation in one is as important as in the other.

If not scaled (and centred) the first gene A will contribute more to the variation than gene B, because the expression variation in absolute numbers is bigger.

Both are used in publications, but I think the scaled and centred is more used, because the first dimension reflects "better" the differences between samples. Of course you can do both and select the one that is better to show what you want.


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

Imagine that I work on an organism without a reference genome, and that the closest reference genome I can get is quite diverged. E.g. ~10% diverged in terms of SNVs when measured with short reads, and also has a lot of structural variants too.

Now imagine I get a 1 million base-pair long-read (e.g. from Nanopore data) for my organism. The question is this:

How can I estimate the proportion of the read that is meaningful sequence vs garbage?

Some things that probably won't work

Most standard approaches won't work here. E.g. I could try mapping the read to the reference, but even if the read was perfectly good I wouldn't expect most of it to map thanks to true structural variations between the read and the reference. The same goes for standard alignment or BLAST.

Some things that might work

The best naive method here seems to be to cut the read up into smaller pieces (either overlapping or not) and use standard approaches to map/align each of these.

So, what have people tried for this? And what tools have you used and why?

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

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


how to search NCBI to find all available reviewed accession numbers for a protein family what is the most optimal keyword search query to input into  NCBI In order to find all the available proteins accession numbers for a specific family?

For example, how would I find all the reviewed accession numbers for all currently available human ALDP (AdrenoLeukoDystrophy Protein) family of proteins that are in NCBI?

Also, is there a specific database within NCBI that is optimized for my task?

Thanks.

A simple approach is to search for the gene name in NCBI Gene portal but this can miss a few cases and often include genes that are unrelated because the gene symbol is used by other unrelated genes as well.

For a cleaner dataset, you can use the NCBI Orthologs as follows. Search for human ALDP in the NCBI homepage with All Databases selected in the dropdown menu as shown below:

After you click the 'Search' button, you will be directed to a result page that looks like this:

Click on the 'Orthologs' button and you will be directed to a page with the list of proteins that are similar to the one you have started from. Note, these may not be orthologs in the truest sense; see this description about how they are calculated.

You can search for the taxonomic group that you are interested in, pick and choose the organisms of your interest and download the table in CSV format. It will have the NCBI GeneID, RNA and protein accessions.


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

cexRow: changes the size of the row label font.

keysize: numeric value indicating the size of the key.

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

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


Transform a R data.frame into a krona plot without krona tools? KronaTools lacks a R pipeline to convert any data.frame into a krona plot. There is a way to go via a phyloseq-class object, but coming from a data.frame I did not find any documentation on how to transform these.

First download any Krona plot (.html) file or make one by using the Excel sheet from the Krona Package. Open in text-editor and remove anything below "datasets" tag and/or the first "node" tag till the end of the document.

Here is the path to go: 1) get your data.frame in hierachical order 2) transform to data.tree object 3) transform to xml file 4) gsub until all the tags match krona annotation

options(scipen=999)
options(stringsAsFactors=FALSE)
library(data.tree)
library(xml2)
library(magrittr)

#get your data into a hieracial format in a data.frame
t <- structure(list(freq = c(34L, 127L, 42L, 123L, 343L, 13L, 104L,
21L, 135L, 91L, 167L, 148L, 90L, 46L, 25L, 79L, 1726L), king.n = c("",
"Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Fungi", "Plantae",
"Plantae", "Plantae", "Plantae", "Plantae", "Plantae", "Plantae",
"Plantae", "Plantae", "Plantae"), class.n = c("", "Dothideomycetes",
"Dothideomycetes", "Dothideomycetes", "Dothideomycetes",     "Dothideomycetes",
"Dothideomycetes", "Liliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida", "Magnoliopsida",
"Magnoliopsida", "Magnoliopsida", "Magnoliopsida"), fam.n = c("",
"", "Cladosporiaceae", "Didymellaceae", "Pleosporaceae", "Pleosporaceae",
"Pleosporaceae", "Poaceae", "", "Asteraceae", "Asteraceae", "Asteraceae",
"Asteraceae", "Fabaceae", "Hydrophyllaceae", "Hypericaceae",
"Malvaceae"), gen.n = c("", "", "Cladosporium", "", "Alternaria",
"Pyrenophora", "Pyrenophora", "Phleum", "", "", "Achillea", "Centaurea",
"Jurinea", "Trifolium", "Phacelia", "Hypericum", "Malva"), spec.n = c("",
"", "cladosporioides", "", "", "", "tritici-repentis", "pratense",
"", "", "", "cyanus", "", "", "", "", "")), .Names = c("freq",
"king.n", "class.n", "fam.n", "gen.n", "spec.n"), row.names = c(1L,
3L, 4L, 5L, 8L, 9L, 11L, 15L, 16L, 17L, 18L, 19L, 21L, 22L, 25L,
26L, 28L), class = "data.frame")

#if any sequences can not be assigned to the highest hieracial level give them a random name like "NA"
t[t[,"king.n"]=="","king.n"] <- "NA"
#define the hieracy (do not include the counting var here), the first space is the root category
t$pathString <- paste(" ",t$king.n,t$class.n,t$fam.n,t$gen.n,t$spec.n,sep = "/")
#transform to a data.tree
t <- as.Node(t)
#define the sequence frequency calculations
t$Do(function(x) x$freq <- Aggregate(x, "freq", sum), traversal = "post-order")
print(t,"freq")
#get all frequencies for each node
repval <- t$Get("freq")

#tranform into a nested list and afterwards into a xml document to get the nesting
t <- t %>% as.list(.,mode="simple",unname=T)
t <- as_xml_document(list(root=t))%>%as.character

#from here it is just a bunch of regex
#split the xml file into one line for each <
temp <- strsplit(t,"\n")[[1]][2] %>% strsplit(.,"(?=<)",perl=T) %>% .[[1]]
temp <- paste0(temp[c(T,F)],temp[c(F,T)]) %>% gsub("[0-9]","",.)
#remove the names nodes and replace by krona annotation
#replace the closing named tags by generic closing "node" tags
temp <- gsub("</(.*)?>","</node>",temp) %>% gsub(">.*$",">",.)
#replace the named opening tags by generic "node" tags with the "name" attribute, replace all (because they are incomplete) numeric values with a placeholder
temp <- gsub("<(?=[^/])(.*)?>",
     '<node name="\\1"><magnitude><val>__\\1</val></magnitude><score><val>1</val></score>',temp,perl=T)
temp %<>% gsub("root"," ",.)
#find the placeholders to insert the frequencies and execute
repos <- grep("__",temp)
for (i in repos)
  temp[i] <- gsub("^(.*)(__.*?)(?=<)(.*)",paste0("\\1",repval[which(repos==i)],"\\3"),temp[i],perl=T)

#add the end of the .html file
temp <- c(temp,"</krona>
","","")  %>% unlist
tout <- paste(temp,collapse="\n")

#append your data to the Krona plot
write(tout,file="Krona1.html",append=TRUE)

Et voilá!


How to tell which channel is the reference channel in two channel (red green) array data? I'm sorry if this is a stupid question - I have not analysed two channel array data before.
I have a GenePix data file, which does not seem to indicate which channel (red or green) is used as reference channel.
If it possible to tell this from the data itself?
Where raw is my data object, and G is a dataframe of the green channel data and R is the data frame of the red channel data, the ranges are as follows:
#range(raw$G)
[1]    45 21939
>
> range(raw$R)
[1]  39 332

No, you can't get this from the data itself. If there is a control/reference sample on one channel, you need someone to annotate which channel it's on.


match two files and check if position of chrom is within exon region I have two files, file1 and file2. I want to match file1 chrom with file2 chrom, and if there is a match then I want to check file2 if the position of the match chrom1 is within the exon region that is in file1. i want to do this with the whole file and print out rows that match and the position is within the exonic regions.

For instance, file1 has a chrom 1 and POS 736689 which match chrom 1 in file2, now I take POS 736689 from file1 and check if it is within the exon region 29553,30563,30975,   30039,30667,31097,... in file2.

If it is within the region it will print out the file1 row that matches the pattern.

The exon regions can be read as:

29553 to 30039
30563 to 30667
30975 to 31097

file1 look like this:

CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO
1   736689  T   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=6
1   758663  CGCC    TTAG    0/0 0/1 0|0 0|1 0|1 0|0 0|1 AC=2;AN=2184
1   762005  G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=13
1   771950  C   T   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=110;AN=2184
1   772172  C   T   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=8
1   790590  C   T   0/1 0/1 0|0 0/1 0/1 0|0 1|1 AN=2184;AC=84

file2 look like this

chrom   Start   End
chr1    29553,30563,30975,  30039,30667,31097,
chr1    53048,54829,    53067,54936,
chr1    69090,  70008,
chr1    160445,161313,  160690,161525,
chr1    324755,327551,  326514,328453,
chr1    367639, 368634,
chr1    459655,461153,461750,   459992,461288,461954,
chr1    523008,529614,  523081,530148,

Despite an awk solution is favored by the OP, here is an R solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward awk solution, at least for me).

res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84


How to run 'join' command for multiple files (of 2 types) in a folder I have two types of files for a gene in a folder

File 1: FOS.tf.txt

ABL2
ACTN4
ADGRE5
ADIPOR2
ADRB2
ERCC4
EZR
FAS
FMO4

File 2: FOS.tt.txt

ABL2
CASP9
CCK
DDIT3
ERCC4
EZR
FAS
FMO4

After running command

join FOS.tf.txt FOS.tt.txt > FOS_final.txt

I got output file: FOS_final.txt

ABL2
ERCC4
EZR
FAS
FMO4

My question is that I have these two type of files (tf and tt) for 605 genes in a folder and I want to run join command together for all my genes and get 605 output files.

If all of your .tf.txt and .tt.txt files are in the same directory, the following command should work.

for tffile in $(ls *.tf.txt)
do
    prefix=$(basename $tffile .tf.txt)
    join ${prefix}.tf.txt ${prefix}.tt.txt > ${prefix}_final.txt
done


fastq file format unknown I have paired-end fastq files some of which seem to be in a weird format (from a collaborator, not a public database). When I cat the file I get what seem to me to be binary syntax but not sure?

However when I cat into the other pair-end the output seems fine

Any help/suggestion/solution would be greatly appreciated.
It does seem to be compressed with a missing suffix:
file pair1.fastq
   pair1.fastq: data.

This indeed looks like some sort of binary file. I have seen bunch of those when people renamed .fastq.gz to .fastq without actually uncompressing them. However, as you say, this is not your case, as file pair1.fastq would identify the file as gzip compressed data.
Looking at possible return values of file comments here, data basically mean the command could not figure out anything about it beside being some sort of binary data file, but I suspect it probably somehow broken file. Honestly, at this point I think you should check the source of the data. If it's from a public database, re-download the data. If it's your own sequencing, contact the sequencing provider and check the data with them. I don't see any way to rescue this file (which does not necerasrily mean there is nothing you can do).