Label protein strucuture residues in Chimera I have generated proteins structure of capsid protein of two strains of viruses. I have modeled them showing Columbic electrostatic potential. I have opened them in two separate windows for comparison. I have following two issues:

1) How can I keep two molecules on the same plane when they are opened in two windows. Is there any function that lets us know which side or which angle the molecule is facing so that I can keep both molecules at same position for comparison.

2) I was able to label residue on the surface. But whenever I move the molecule some of the residue start moving and some remain static. I want to show that there is difference is electrostatic potential because of difference in the residue at some particular site. I have not been able to do that because those residue keep moving haphazardly?

It would be great if someone could help me with this.

You cant do this in Chimera, but in pymol you can write a script to do this.

In terms of the view, you can use get_view on the command line. This returns a section that you can copy and paste in - which will give you the same view every time. I recommend labelling the residue you want, then finding the position you want to view it from. Then using get_view to get that view, which you can use every time (just copy it into a text file) for a reference

This does presume your different molecules are aligned, but you can do that in pymol as well to make sure that when you use the output of get_view it is consistent.


How to combine multiple kernels of large sample datasets? I have multiple large sample datasets in matrix format (each has 15000 rows and 5-50 columns) corresponding to different experiments. Each matrix contains the same number of samples(rows) but the variables(columns) are not the same. My objective is to cluster the samples on the basis of all the experiments.

I tried to use Unsupervised multiple kernel learning (UMKL) to integrate the datasets followed by kernelPCA using "mixkernel" package in R (https://www.ncbi.nlm.nih.gov/pubmed/29077792). The UMKL step calculates kernels for each dataset and combines the kernels using 3 different approaches: 1) calculating a consensus kernel from multiple kernels 2) calculating a sparse kernel preserving the original topology of the data 3) calculating a full kernel preserving the original topology of the data

The kernel calculation step was fine but the kernel integration step (all 3 approaches) runs very long and my computer hangs.

Is there any way to handle this problem? More specifically, is there any way to handle multiple kernel integration for large sample datasets?

Any suggestion alternative to using kernel methods will also work.

Welcome to Bioinformatics Stackexchange @SD1024. A recently published algorithm known as MOFA (Multi-Omics Factor Analysis, paper, github) is generating a lot of interest, and is designed answer exactly the sort of biological question you are describing. It claims to extract axes from multiple matrices with overlapping samples, but not necessarily overlapping features and handle missing data elegantly. I havn't tried it myself you, but I'm quite excited to do so.


How to understand and analyse RNA-seq data (for a beginner)? I am trying to understand expression of a certain protein across Pseudomonas species. I downloaded an SRA file from NCBI and converted it to a fastaq file. I am not able to understand how to interpret and analyse this file. I want to get the mRNA sequence and see which genes are being expressed and find the mRNA of interest.

How do I proceed with getting the mRNA sequences and then for mapping them to genome?

SRA files are raw sequencing data that have been converted to the NCBI storage format. These can be converted to fastq files, e.g. with fastq-dump from the SRAtoolkit.

Beyond that I suggest you start with the Hitchhiker's GUide to RNA-seq and the DESeq2 workflow to get a background. Be sure to extensively browse the web for vignettes, guides and tutorials. YOu will find helpful resources to get started.


I am getting this error when running a test data set for Lyve-SET on RHEL 7. Is this an MPileup error or something else? This is the error that I am getting:

launch_varscan.pl: Creating a pileup /home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.mpileup
removed ‘/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.1239.mpileup’
launch_varscan.pl: Running mpileup:
  echo "NC_001416.1
" | xargs -P 1 -n 1 -I {} sh -c 'echo "MPileup on {}" >&2; samtools mpileup -f lambda/reference/reference.fasta
-Q 0 -B  --region "{}" /home/user/lyve-SET-1.1.4f/lambda/bam/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam >
/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.$$.mpileup'
MPileup on NC_001416.1
mpileup: unrecognized option '--region'
open: No such file or directory
[mpileup] failed to open NC_001416.1: No such file or directory
launch_varscan.pl: main::mpileup: ERROR with xargs and samtools mpileup
Stopped at /home/arahat/lyve-SET-1.1.4f/scripts/launch_varscan.pl line 106.
QSUB ERROR
31488
launch_set.pl: Schedule::SGELK::command: ERROR with command: Inappropriate ioctl for device
  /usr/bin/perl /home/user/lyve-SET-1.1.4f/lambda/log/SGELK/qsub.2705.pl
Stopped at /home/user/lyve-SET-1.1.4f/scripts/../lib/Schedule/SGELK.pm line 661.
set_test.pl: ERROR: died while trying to run COMMAND
  launch_set.pl lambda --numcpus 1 --numnodes 20
set_test.pl: See if you have the correct software installed for this command.
6400 at /home/arahat/lyve-SET-1.1.4f/scripts/set_test.pl line 94.
make: *** [test] Error 25

I think the main problem is the "mpileup: unrecognized option --region", which seems to be blank. In the beginning of the block, it says creating pileup and then it immediately removes the file, probably to create the pipeline file.

I cannot change the code for Lyve-SET pipelining or create shell scripts for the pipeline because I am only configuring my server for use with the genome sequencing pipeline. Would this be a problem with the data set or some tool that I need to install? Any insight would be appreciated.

The samtools v0.1.19 manual confirms that there is no --region option.  In that version there is only the -r argument to pass in the targeted region.  If you update to a newer version of samtools this should execute correctly


BiomaRt error: Error in martCheck I am experiencing an error while using getSequence() in biomaRt.

The code is given below:

library("biomaRt")
ensemblp <- useMart(biomart="plants_mart",host="plants.ensembl.org")
osj = useDataset("osativa_eg_gene", mart = ensemblp)
genes <- read.csv("osj.csv")
osj_id <- genes$ensembl_transcript_id

os_genes <- getSequence(id = osj_id, type = "ensembl_transcript_id", seqType = "peptide", mart = ensemblp)

The error message I am getting is :

Error in martCheck(mart, c("ensembl", "ENSEMBL_MART_ENSEMBL")) :
This function only works when used with the ensembl BioMart.This
function only works when used with the ENSEMBL_MART_ENSEMBL BioMart.

Did I make any mistake in code? Please help me fix the error. I tried asking this in Bioconductor support, but no one responded yet.

getSequence has only been enabled for the main Ensembl (vertebrates) biomaRt. This is all at the Bioconductor end. You'll need to use a getBM instead, for example:

os_genes <- getBM(attributes =c("ensembl_transcript_id", "peptide"), filters = "ensembl_transcript_id", values = osj_id, mart = osj)


How to confirm exon shuffling in a gene? Note: this question has also been asked on reddit

I'm trying to confirm that the sequence of a novel gene is derived by exon shuffling between several different genes. I have the promoter sequence, gene sequence, and mRNA (with defined exon/intro boundaries). I've tried performing several database searches with the sequences, but each search produces a set of different unrelated hits. How would one go about using this info to confirm the exon shuffling hypothesis? Any tips are appreciated. Thanks!

These suggestions are based on your comment that you have the gene sequence.

For a statistical approach, you could formulate a null model and show that the sequence similarity you observe is much higher than expected under the null model, but reasonable under your exon shuffling model. Ideally, you would show this individually for each exon, to convince the reader about the exons' disparate origins. I'll propose a reasonable null model, though suggestions are gladly accepted here. Pick a random gene, assume it's built using exon shuffling, and search the rest of the transcriptome for sequences that best substantiate that argument.

For an evolutionary approach, you could search for orthologs (homologues? I never remember the difference) of your gene in other species. Hopefully, the results would show that your gene formed more recently than the exons it's built out of. If your gene actually has some other origin, maybe you'll find out it's older than some of its putative components. This can't exactly "confirm" your hypothesis, but it will test your hypothesis and just be generally a cool thing to try.


About CAFA experiment I have a question about CAFA1 and CAFA2 experiment.

CAFA is an experiment that compare computational methods predicting protein functions. In the paper discussing CAFA2, there is an image comparing 2 methods (CAFA1 versus CAFA2):

How do I interpret the numbers and colors (dark and light)? What do they mean?

Note that in the article published that figure is different from the preprint!
See the figure:

As they aren't substantially different in display other than showing the two set of predictions that are comparable between CAFA1 and CAFA2 I explain them below:

The colors are the margins of $F_{max}$ they refer to the equation found in the section "Comparison between CAFA1 and CAFA2 methods" under methods:

$$\delta(m_2, m_1) = \frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_2)-\frac{1}{B}\sum_{b=1}^{B}F_{max}^{b}(m_1)$$

Where B is the 10000 bootstraps they did and $m_1$ and $m_2$ are the CAFA1 and CAFA2 methods.

Where $F_{max}$ is explained in the article (and in the preprint) as:

  A perfect predictor would be characterized with $F_{max} = 1 $

So, the higher the (positive) difference is, the better precision and recall (See wikipedia) are for those methods. Note that the scale is from -0.1 to 0.1 not from -1 to 1 so the differences are quite small (One could argue that the log10 would be better to see those small differences)

One can appreciate that in the bar plots at the right, where the differences between each dataset are small in the B section, so the color is white, and in the C panels they have brighter color.

The numbers in the boxes are the percentage that a CAFA2 method did better than a method in CAFA1 for another dataset in the "intersection of the two GO snapshots" as noted in the methods section. (But it is also not clear to me)


alternatives to MEDIPS to analyse MeDIP datasets MEDIPS is an established tool with functions for the quality control and analysis of data derived from immunoprecipitation (IP)-seq samples, like Methylation IP sequencing datasets.

I would like to know if there are any other tools I should consider nowadays that could do a better job and hence replace MEDIPS for the analysis of MeDIP datasets.

I guess QSEA should be an obvious answer:

  QSEA (quantitative sequencing enrichment analysis) was developed as
  the successor of the MEDIPS package for analyzing data derived from
  methylated DNA immunoprecipitation (MeDIP) experiments followed by
  sequencing (MeDIP-seq).

See http://bioconductor.org/packages/release/bioc/vignettes/qsea/inst/doc/qsea_tutorial.html


How to remove zero value of gene on FeatureScatter plot using Seurat? Nowadays, I am trying to calculate Pearson correlation values between two genes of my interest from single-cell RNA-data (features.tsv, barcodes.tsv, and matrix.mtx files) which are obtained from the 10xgenomics cellranger. I am interested in tumor-infiltrating Tregs. I downloaded fastq files from a published paper that describes tumor-infiltrating Tregs with IL10+IL35+ (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126184), and ran "cellranger count".

Yesterday, I learned how to get a Pearson correlation value from a kind guy. I made a graph below.

In general, Tregs are defined by CD4 and Foxp3 expression. As a trial, I want to find the correlation between two genes. The correlation value was 0.05, which might have been due to the presence of Tregs that have zero value of either CD4 or Foxp3. So, my question is that can we remove cells with 0 value of genes to get more high correlation value? My data looks like this.

> Treg.data[c("Cd4", "Foxp3"), 1:30]

------------------------------------- # results

2 x 30 sparse Matrix of class "dgCMatrix"
   [[ suppressing 30 column names ‘AAACCTGAGAAACCAT-1’, ‘AAACCTGAGAAACCGC-1’, ‘AAACCTGAGAAACCTA-1’ ... ]]

Cd4   . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Foxp3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

------------------------------------- # results

Each dot represents 0 value.

Meanwhile, I performed to standard pre-processing workflow, as suggested in Seural Tutorial. For example,

> Treg[["percent.mt"]] <- PercentageFeatureSet(Treg, pattern = "^MT-")

> Treg <- subset(Treg, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

Can you help me? Thank you very much.

It will not improve your correlation most likely, unless the zeros are not missing at random:

library(Seurat)
data <- pbmc_small
data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^MT-")
data <- NormalizeData(data)
genes <-c("HLA-DRA","LYZ")
countdata = pbmc_small[["RNA"]]@data
countdata <- t(as.matrix(countdata)[genes,])
plot(countdata)

Now you need to subset the matrix, and remove the rows where at least one of the entries is 0, you can use rowSums(..==0)==0, meaning no zeros in the rows:

plot(countdata[rowSums(countdata==0)==0,])

cor(countdata[rowSums(countdata==0)==0,])
          HLA-DRA       LYZ
HLA-DRA 1.0000000 0.2338371
LYZ     0.2338371 1.0000000


What does that mean if there is no bifurcations in my time series single cell data? Sorry,

I have analysed my time series single cell data with Scuba algorithm but that says that there is no bifurcations in my data and outputted a single chain of 4 states with no branching like below figure

Sorry, I googled but I am not quite sure what this means

Any suggestion please?

Thanks a lot

The most obvious interpretation is that the four cell lines are related in a linear pathway rather than a branching one.


How to load and split my two individual datasets when integrating datasets in Seurat? I want to do comparative analyses of my two individual single cell RNA-Seq datasets in Seurat. And I load and list my data as follows:

Control.data <- Read10X("~/Documents/...  Control_filtered_feature_bc_matrix")

Control <- CreateSeuratObject(counts  = Control.data, project = "Control_summary", min.cells = 3, min.features = 200)

Mut.data <- Read10X("~/Documents/... /Mut_filtered_feature_bc_matrix")

Mut <- CreateSeuratObject(counts  = Mut.data, project = "Mut_summary", min.cells = 3, min.features = 200)

Object<- list(Control,Mut)

names(Object) <- colnames(Object)

Then the protocol https://satijalab.org/seurat/v3.1/immune_alignment.html said, I need to use SplitObject function, but what do I split by?

Object_list <- SplitObject(Object, split.by = "?")

Seurat is expecting individual datasets to be normalized separately prior to data integration. In that respect lists of objects corresponding to different datasets are handy to manipulate each object/dataset individually.

In your case you do not need to do any splitting, you have an object corresponding to each of your datasets.

If you would have started with a single count matrix of all datasets then you would need to split your data based on batches/experimental condition/...


Infer the new (Ban et. al) ribosomal nomenlature (ex. uL53 ) from the ribosomal protein's sequence I'm trying to programmatically construct a name for each protein subchain in any ribosome from Uniprot in accordance with Ban et. al's 2014 proposal (excerpt given below) using PDB's and Uniprot's APIs.
This is to be a part of a bigger ribosomal analysis suite released next year and we would like to standardize the names as much as possible.
Ban et. al's proposal: [https://bangroup.ethz.ch/research/nomenclature-of-ribosomal-proteins.html]

I'm however not remotely a biologist (I work in compilers) and would really appreciate some clarification on the way nomenclature works across domains of life and how homology plays into it.
This is my naive approach thus far:

I get a ribosome from PDB let's say (3J9M or 5MYJ or any other for that matter) and split it into protein subchains. So on the order of 20-80 proteins.
For each protein subchain, I can get all the information in the world on it from Uniprot: names, sequence, organisms containing it, publications, synonymous names, you name it...
Right now (again, naively) for those proteins that don't already have a new name assigned to them, I take name-synonyms that PDB has plenty of (from various publications, I assume) and scan them with regexes (ex./[LS]\d{1,2}/g) for things like "L8" or "S15".
Thus, for each subchain, I end up with something like "L15" or "S28" in multiple copies (depending on the abundance of synonymous entries in PDB).
The rut i'm in right now is this: having something like "L15" to match it to a new "Ban-nomenclature" name: they provide the nomenclature look-up tables in terms of these three "historical" naming conventions(i.e "human", "yeast", "bacteria") and a taxonomic range while all the I have is the ribosomal structure(ex `` from which I got the subchain. I'm sure there's plenty of information that I can get

Could somebody elucidate to me which assumptions I can and cannot make in this field and how exactly to use the taxonomic range?

Say, (i) can I assume that all archean names are the same as bacterial homologs?
If I draw a eukaryotic cytoplasmic ribosome, do I match its subchains against yeast column or human column (since both are eukaryotic)?

Somebody suggested that I start grouping these subchains based on sequence, not on a domain and I do have access to each protein's sequence, but I am not sure how to go about it. I was thinking to get each of Ban's new names' UniRef90 cluster, let's say, and then start checking which cluster a given subchain belongs to.

But I'm obviously guessing at this point. Any advice or pointers are very appreciated.

This is a copy-paste nearly verbatim of comments so that the question has an answer

What you are trying to do is find what genes are in what homology cluster. This is common problem and there are many solution each with some issues.

Uniprot90 is indeed a cluster of homologues but it is too limited. Whereas you require clusters that span all the domains (universal).

COG

The nearly original grouping, is COG (cluster of orthologous genes). It is really handy for this but aren't really available from Uniprot and NCBI (its curator) has been phasing COGs out. In Uniprot records, they are patchily present as its derivative eggNOG (which has varying degrees of "scope", LUCA being the universal one):

<dbReference type="eggNOG" id="COG0451">
<property type="taxonomic scope" value="LUCA"/>
</dbReference>

Pfam

Another option that is well mapped and always present when possible in Uniprot entries is Pfam ids. These are domain family annotations and are also good option. It is a domain fold group. So domains within protein with homologous folds are grouped into one. If a protein has two known domains, it will have two Pfam entries —repeat proteins have loads, but most protein have one domain. In Uniprot XML you will find these are:

<dbReference type="Pfam" id="PF01370">

Similarly there is InterPro, which is more narrow than Pfam, but generally overlaps.

DIY

Then there's always the option of making your own psi-Blast motif on the NCBI website (setting high number of hits, setting the database to refseq or even PDB and doing a few iterations), saving the motif pattern (PSSM) and using local blast with it.

Caveats

In general you may however encounter these pitfalls:

you may have two protein clusters that share a single id
you have gene fusions (one protein, 2 domains of interest). In Pfam the protein will have two ids, with COGs or similar, you'd have only an anamalous length to go by.
you may have wished to find analogues (i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural). No homology clustering scheme can find these.  Say a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).
you have a very small protein which is poorly conserved


Any fast options to query large VCF bed intervals? I'm doing some analysis and I need to subset a large VCF file (~8GB gziped) given a bed interval and identify within a subset of rsid.

Unfortunately, both my normal choices to do this analysis (snpSift and bedtools) are taking way to long or failing due to memory issues in my local computer and a remote server.

Do you guys know any other options or suggestions to speed up this process?

Follow the commands I use:

bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

or

gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

The bedtools command usually fail due to unknown reason and SnpSift runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.

this:

bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.


Difference between CPM and TPM and which one for downstream analysis? What the difference between TPM and CPM when dealing with RNA seq data?

What metrics would you use if you have to perform some down stream analysis other than Differential expression for eg.

Clustering analysis using Hclust function and then plotting heat map to find differences in terms of expression levels, correlation and pca

Is it wrong to use TPM for such analysis, if yes then when does one use TPM versus CPM.

You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-normalized values of the other genes).
If one has to choose between those two choices one typically chooses TPM for most things, since generally the length normalization is handy. Realistically, you probably want log(TPM) since otherwise noise in your most highly expressed genes dominates over small expression signals.


Is there public RESTful api for Gnomad? I currently find Harvard's RESTful API for ExAC extremely useful and I was hoping that a similar resource is available for Gnomad?

Does anyone know of a public access API for Gnomad or possibly any plans to integrate Gnomad into the Harvard API?


Which are the use cases for the methods for DE in Seurat In Seurat we can specify multiple methods for finding DE genes.

I am wondering when should we use which one, the use cases. As I understand, whole Seurat package is used for single cell sequencing analysis. I ran different methods for my dataset that was produced from brain cells, around 800 cells in each sample, by single cell sequencing. Now I am confused which methods are the most appropriate to use, which results are the most reliable. For single cell sequencing when would you prefer to use one method over the other?

You can take a look at the recently published article: Bias, robustness and scalability in single-cell differential expression analysis.

  We evaluated 36 approaches using experimental and synthetic data and found considerable differences in the number and characteristics of the genes that are called differentially expressed. Prefiltering of lowly expressed genes has important effects, particularly for some of the methods developed for bulk RNA-seq data analysis. However, we found that bulk RNA-seq analysis methods do not generally perform worse than those developed specifically for scRNA-seq.

Then you can decide which method fits better to your type of scRNA-seq protocol. Since Seurat offers DESeq2 could be a good idea to use it, but this depends a lot on your data.


Parse BAM by insertion size and get genomic coordinates I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, file.bam.

(In terms of genomic coordinates, I just mean I would like to have a rough idea where to look at the BAM in IGV, chromosome and position.)

Based on the BAM CIGAR, I should be able to grep the BAM based on some cut-off of "large inserted" bases.

It's unclear to me how to best do this via samtools; one could parse the cigar string via samtools for insertions, but I don't understand how to translate this into the genomic position in the alignment.

This goes with various samtools wrappers, e.g. pysam in python:

import pysam

bamfile = pysam.Samfile(bamFile, "rb");

for read in bamfile:
    if not read.is_unmapped:   ## only get mapped reads
        cigarLine = read.cigar
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1: ## this means it's an insertion...
                  ## somehow filter only large insertions, and then get coordinate information

I'm not sure how to do this in a non-clumsy manner. Any help appreciated.

You basically have the right framework already in your pysam code, which I took the liberty of tidying up a bit. You just need to add a bit:

for read in bamfile:
    if not read.is_unmapped:
        cigarLine = read.cigar
        interesting = False
        for cigarType, cigarLength in cigarLine:
            if cigarType == 1:
                if cigarLength > 100:  ## Arbitrary value
                    interesting = True
                    break
        if interesting:
            print("{}:{}-{}".format(read.reference_name, read.reference_start, read.reference_end))

Unless you have long reads, this will be close enough for IGV.


Assign cell types to groups of cells based on their gene expression profiles I have large filtered, normalized dataset of scRNA-seq data of C.Elegans species. Rows are genes (10 000), columns are cells (66 000). Let's say that I got 40 different groups of the cells based on their expression profiles, how could I now assign cell types to these groups? I would expect to get some reference dataset for C.Elegans and compare it with what I got, but I honestly have no idea of where to get such a dataset, and whether there are any R packages that do that. Any advice would be greatly appreciated.

I think the analysis will involve two steps:

1) Compile a reference list of cell types, with a list of genes expressed in each of them, based on the literature. C. elegans has only ~1000 cells, so I'm sure it's been well researched. Mark which genes are just expressed in that particular cell type, and which are uniquely expressed, this latter will be more useful markers.

Resources to do that:

http://openworm.org/science.html#other
http://www.wormatlas.org/
http://www.wormbase.org
Pubmed.
Your lab's knowledge.

This will involve mostly manual work, but maybe you can find some list.

For example, platelets express VWF. And activated platelets express SELP.

2) Use one of the single-cell RNA-seq packages (e.g. Seurat) to cluster cells and to find the "marker" (DE) genes for each cluster (vs. all other clusters).

For each cluster, look up the marker genes in the reference list and infer the cell type.

You may need to infer the cell type from the functions and descriptions of the genes.
Also see Retrieve detailed gene descriptions

Keep in mind that you will have many confusing doublets in the dataset expressing characteristics (genes) of two cells.

You may also need to re-do t-SNE and clustering on a subset of cells, e.g. blood/immune cells, in an iterative, hierarchical way.

Do you have information on which cells are from the same worm? Once you have a reference list, you can try and identify cells of just one worm.


Spearman correlation between two genes I have FPKM data for three genes like below:

     Samples      gene1 gene2 gene3
TCGA-2Y-A9GS-01A   0.9   7.4  1.0
TCGA-2Y-A9GT-01A   0.8   1.0  0.3
TCGA-2Y-A9GU-01A   0.6   2.0  0.2
TCGA-2Y-A9GV-01A   1.2   0.5  0.1
TCGA-2Y-A9GW-01A   3.8   2.1  0.4
TCGA-2Y-A9GX-01A   2.3   2.0  1.5

I used cor.test

cor.test(~ gene1 + gene2, data = df2, method="spearman", continuity=FALSE, conf.level=0.95)

    Spearman's rank correlation rho

data:  gene1 and gene2
S = 5686100, p-value = 6.083e-09
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho
0.2984045

I have a warning message which I didn't see before.

Warning message:
In cor.test.default(x = c(0.9, 0.8, 0.6, 1.2, 3.8, 2.3, 3.8, 0.4,  :
  Cannot compute exact p-value with ties

Do I need to care about this warning message? Is it good using FPKM data for correlation?

For plotting I used ggscatter.

ggscatter(data, x = "gene1", y = "gene2",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "spearman",
          xlab = "gene1", ylab = "gene2")

Is this fine or I need to use any log for the scatter plot?

Since Spearman is a rank-based test, it relies on you being able to accurately decide on the ranking of your observations by some metric (usually the magnitude of the numbers). If two observations have identical values (are tied), then they cannot be definitively ranked. Since the ranks are not unique, exact p-values cannot be determined.

e.g. in your data, for gene 2:

                  gene2 rank
TCGA-2Y-A9GS-01A    7.4    1
TCGA-2Y-A9GT-01A    1.0    5
TCGA-2Y-A9GU-01A    2.0   3=
TCGA-2Y-A9GV-01A    0.5    6
TCGA-2Y-A9GW-01A    2.1    2
TCGA-2Y-A9GX-01A    2.0   3=

By way of example, it's possible to do the test with either TCGA-2Y-A9GU-01A or TCGA-2Y-A9GX-01A ranked number 3 - meaning a definitive p-value cannot be determined.

If there are only a few ties, you probably don't need to worry too much. If your FPKMs (or better, TPMs) were calculated with greater precision, then ties would be less likely. Is there a reason you only have these numbers to 1 decimal place?

Log or linear scale doesn't make a huge difference (and neither is 'wrong'). If you choose to plot log counts, you need to consider what happens if you have zeroes in your data.


Automating download of results of Drug and Drug Target Mapping The drug and drug target mapping tool shows the top 3 PDB IDs for each drug based on the drug target sequence similarity search and ligand id search. Is there a way to get all the PDB ids of the drug-target PDB files and not only the top 3?

I know you can get this information after clicking the "find" button under the last column of the output table (the column named "All PDB"), but I have a query of several hundreds of drugs and I'm not willing to do it manually for each result.

Simpler option would be to not use the query interface  and directly parse the CSV file found on that page.


List of global ethnicities Looking at the 1000 genomes project I see the following list of population groups and super-groups:
https://www.internationalgenome.org/category/population/

However, is there something like a standard list of global ethnicity? Like an ontology of population diversity or similar?

I'm trying to tie genome data to 'ancestry', and to do that I'd like to survey a set of people on their ancestral background (as far back as grandparents ideally).

How to capture diversity in terms of what people 'self report' in their backgrounds? e.g. middle / far east, north / south India, etc.

It's not clear exactly what your goal is. It might help to have an idea of how you are trying to use this data; and it's probably not a bioinformatics question, as ethnicities are social/political rather than biological constructs.

However, there exist some controlled vocabularies used by LoC for ethnic groups; that looks a bit old so it might not be the best source.

I am by no means an anthropologist but it does look like there are anthropological databases devoted to e.g. ethnography that may have complementary information.

Wikipedia has a simple list of ethnicities that notes, instructively, "There has been constant debate over the classification of ethnic groups."

I think that you are unlikely to find a single authoritative answer that is safe to use as a source of truth. I would normally suggest to try another stackexchange site, though it looks like there is none for sociology or anthropology. There is no anthropology tag available on this site, but you might try posting in e.g. history.stackexchange.com and adding an anthro tag.


Waterman-Smith-Beyer Sequence Alignment I am having trouble understanding the affine gap penalty in the following example -

I am not sure where the 3 and 4 come from or the 4 and 5 in cells 1,3 and 2,4. I'm not understanding how the affine gap penalty works here. I understand how the linear gap penalty worked in the Needlemen_Wunsch but not this.

The linear gap penalty has a cost which is directly proportional to the length of the gap: the gap costs 1 unit per consecutive missing base.
The simplest case is assigning each unit (= each missing base) a penalty of 1:

  g(k) = k

Or more generally, if you want to increase the penalty per base to 2, 3 or n:

  g(k) = n * k
  with k = length of the gap, and n = cost of 1 missing base

The affine gap penalty has a cost that follows -as its name conveniently states- an affine function.
The idea behind it, is that any gap in an alignment has a cost with two components:

the "gap opening" penalty: it is a cost payed only once per gap, basically just for existing.
the "gap extension" penalty: equivalent to the linear gap penalty, it costs 1 unit for each consecutive missing base.

The general formula for the cost is then:

  g(k) = a + n*k
  with a= opening penalty, n= extension penalty, k= length of the gap

Compared to the linear gap penalty, this will give a small bonus to longer consecutive gaps compared to skipping one nucleotide, matching one, skipping one again etc (because in the case of longer gaps, you only pay the opening penalty once, so the global cost per base decreases).

Cells on your diagram:
From what I understand, you're working with an affine penalty of g(k) = 1 + 1*k. For the [3,4] you're "coming" from the score of 2 on the left (cell 1,2). You either:

match the A (row 1) with G (column 3) and get a penalty of 1 for substitution -> total of previous 2 + new 1 = 3
open a gap to skip the G (column 3) and get a penalty of 2 (1 for opening + 1 because gap length of 1) -> total of previous 2 + new 2 = 4.

So you choose the lesser penalty and go with 3, which is the number selected in the bigger box.
The explanation is the exact same one for the other cell you mentioned.


replacing SNPs with missing calls with a specific string I have a big file containing 27 columns and nearly 6 million rows. The following is a little example of my file

head data
0.65   0.722222   1.0      0.75     0
0.35   0.277778   0.0      0.25     0
0      0.666667   0.75     0.5      0.5625
0      0.333333   0.25     0.5      0.4375

Rows are SNPs and columns are individuals. and I have "2 rows per each SNP" (one for reference allele and the other for alternative allele). In the example above I showed data for 2 SNPs (rows 1 and 2 correspond to SNP1 and rows 3 and 4 correspond to SNP2). I want to check if for each SNP, both ref and alt allele frequencies are 0, replace them with 9. this is my desired output:

head desired
0.65   0.722222   1.0      0.75     9
0.35   0.277778   0.0      0.25     9
9      0.666667   0.75     0.5      0.5625
9      0.333333   0.25     0.5      0.4375

A possible approach in plain python:

#!/usr/bin/env python3

import sys

def both_zero(freq1, freq2):
    return freq1 == "0" and freq2 == "0"

with open(sys.argv[1], "r") as snp_file:
    for line1 in snp_file:
        # Consume one extra line to get the frequencies
        # for the alternative allele
        line2 = snp_file.readline()
        # "zip" the frequencies in pairs
        freq_pairs = zip(line1.strip().split(), line2.strip().split())
        recoded_pairs = [
            ("9", "9") if both_zero(freq1, freq2) else (freq1, freq2)
            for (freq1, freq2) in freq_pairs]
        # "unzip" the recoded frequencies into separate objects
        # (see https://stackoverflow.com/a/19343/1878788)
        (recoded1, recoded2) = zip(*recoded_pairs)
        print(*recoded1, sep="\t")
        print(*recoded2, sep="\t")

Calling it with your example snp frequencies file as command-line argument:

$ ./recode_absent.py snps.txt
0.65    0.722222        1.0     0.75    9
0.35    0.277778        0.0     0.25    9
9       0.666667        0.75    0.5     0.5625
9       0.333333        0.25    0.5     0.4375


Violinplot of gene expression I have plotted the log normalized expression of two genes by violonplot for 4 clusters. I have links to my pictures and Seurat object too. My problem is this; in violin plot I can not see the mean or any centennial tendencies so that I don't know if two genes is expressing higher or lower in contrast to each other in each cluster.

violonplot vlnPlot(seuset,c(""))

Actually when I am thinking I see I need a violin or box plot showing  the  ratio of these genes to each other in each cluster and how nice to have mean also in the plot

Somebody knows how I can do that?

You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.

(1) First, notice that vlnPlot() is deprecated. Use VlnPlot().

You can try using the parameter do.sort=T:

VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)

Alternatively, you can return the ggplot2 object and then plot the means etc:

p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p

(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.

In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.

Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group.
(a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.

Here are some functions for retrieving and plotting data from the object:

ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=1)

stripchart(datatoplot2,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=2)


Occupancy of TFs with the target genes The occupancy of SMARCD3 in the target genes listed below. I want to see average, normalized ChIP-seq signal at the promoter proximal region (1000bp upstream and downstream of the TSS).

I have 4 different experimental conditions (overlayed in one plot, per target gene using different color codes) where we want to visualize the changes in occupancy of GATA1 under all the conditions at individual gene promoters. The following are our target genes

ABCB1
ABCC1
ABCC2
SNAI2

So how GATA1 occupancy is changed across all these target genes across 4 different treatment condition.

So far what i have done

Aligned the samples
Done peak calling
Annotated the peak

I did steps 2 and 3 using homer.

Now after all this how do I find the occupancy of a TF with its target genes? Do I get the information from my peaks annotated output or from the peak files?

I would ignore peak calling for this and instead compute enrichment of ChIP/input for the genome (e.g., with deepTools or presumably homer) and then plot it for the genes of interest individually (e.g., using IGV or pyGenomeTracks) or as a group (e.g., with computeMatrix). If the peaks are obvious and you trust your peak calling then sure you can just use that instead, but given that that's prone to error it's usually better to simply use normalized signal.

Regarding actually looking at occupancy, I don't think that's useful. Occupancy isn't binary. If you're instead wanting to associate peaks with genes then using the nearest gene is common (it's not great, but it's not terrible either).


SRA Toolkit and lebanese data I am trying to extract data from this:
https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR6245218

I installed SRA Toolkit, downloaded the SRR6245218 file and executed this:

fastq-dump -X 5 SRR6245218 --fasta -O ./

The result is a fasta file 900 bytes long, which does not contain the mitochondrial data that I need.
What am I doing wrong?

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


How to obtain amino acid context at a given position I'm trying to find the amino acid context around a particular position in a protein (5 aa amino acid window on each side).

I used to have some code that works and then biomaRt stopped working. I think some of these issues would be fixed in later versions of biomaRt, but I would have to install a new version of R, which would require me to install a new version of my OS.

require('biomaRt')
require('Biostrings')
hs_ensembl_mart =
useMart(biomart="ENSEMBL_MART_ENSEMBL",dataset="hsapiens_gene_ensembl")

gene_names=c('KRAS')
aa_affected_loc=12

biomart_query_results=getBM(attributes=c('hgnc_symbol', 'entrezgene','coding'),
  filters = 'hgnc_symbol',
  values = gene_names,
  mart = hs_ensembl_mart)

preceding_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc-5,aa_affected_loc-1)

trailing_peptide= substr(as.character(translate(DNAString(biomart_query_results[1,'coding']))),aa_affected_loc+1,aa_affected_loc+5)

Does anyone have a workaround for this?

It turned out that I needed to reinstall OSX to upgrade my R version and biomaRt install.


Repetition-free Longest Common Subsequence as similarity measure I'm starting a thesis in Computer Science with topic the RFLCS (Repetition-free Longest Common Subsequence). I know very little about Bioinformatics problems and Genomics. I would like to understand the motivations behind this problem.
More precisely, I'd like to understand why one need a similarity measure between genomes, and how the RFLCS is an useful measure.

[Adi et al.] introduce the problem in this article.

A brief literature review can give you the answers.

This preprint gives the definition of RFLCS more clearly:

  A repetition free Longest Common Subsequence (LCS) of two sequences x and y is an LCS of x and y where each symbol may appear at most once.

Adi et al proposed to solve the genome rearrangement problem raised in this paper. If you want to know more, follow the references and citations of these papers and read them.

Basically, Adi et al think the length of RFLCS measures the distance between two genomes. The longer the more similar. I am not sure if this is a good measurement, though.


Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`? I am constructing a bit of software which pipes the outputs of the bam file via samtools view into a script for parsing. My goal is to (somehow) make this process more efficient, and faster than samtools view.

I am only using 3-4 fields in the bam. Naturally, if I only wanted those fields, I could do something like samtools view file1.bam | cut -f ##.

However, I'm trying to be faster than this approach. More specifically:

(1) Given the bam file format, is it theoretically possible to output several requested SAM fields faster than samtools view|cut -f?

(2) Is there software available to achieve (1), or would one need to hack samtools?

The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power.

As this question suggests, if only a specific field from alignments is needed (or fields), then it will probably be better to extract just those fields and do any necessary conversion only on those fields. While this can be done by writing a BAM parser from scratch, many developers have already written software libraries to process BAM files in this way.

Devon Ryan has suggested htslib, which is a C library written by a group that includes the people who wrote the SAM/BAM/CRAM file format specifications. There's also pysam, which is a python wrapper around htslib.

The particular tool that is used will depend on your familiarity with programming and the specific thing that you want to do. If you just want to "output a BAM file", then cat is the quickest:

cat file1.bam

... but you probably don't want to do that, because it seems like you want to process a text-based representation with a script. Because you haven't specified which fields you are interested in, it's not possible to suggest the best thing to use. In the end, I expect that an "efficient" solution to your problem would involve htslib in some form. This is not really hacking samtools, it's using the backend of samtools to process BAM data.

However, it's useful to also think about the time cost of coding. What application do you have which means that the text-based processing method is not sufficient? It takes time to write code, and a lot more time to debug that code to make sure it's doing the right thing in all situations. If this is for a one-off thing, then samtools view output fed into your script may be the quickest solution.


How to construct the phylogenetic tree of the whole genus? I am confused about how to construct the phylogenetic tree of the streptomyces genus. The taxonomy of the streptomyces genus is mostly based on the MLST method, but now I need to use the gene data from NCBI to construct the phylogenetic tree of the streptomyces genus. Should I use the whole gene of the MLST's housekeeping gene to construct the tree? Or I should just use the loci to build the tree? Or the other method? I am a junior in university, have little knowledge of the filed., but my tutor asked me to do that. plus: Sorry for my awkward expression~

The MLST genes are not used to contruct a tree, because each pattern of alleles determines a predetermined genotype, e.g.  Staphlococcus aureus ST 22 is the MLST type concerning the MRSA epidemic. You could build a tree, but its only purpose is to visualise the ST, so its a bit pointless but okay at least you will not mistype something.

Performing a whole genome tree on Strep for the genus is complicated but doable and there are missing genes and stuff and if you just want typing then its not needed.

So if typing is your goal I would submit your MLST loci to the Steptomyces MLST database,

https://pubmlst.org/bigsdb?db=pubmlst_streptomyces_seqdef


How does picard's MarkDuplicate handle unmapped reads? Our BAM files are created according to a "lossless" alignment procedure [1] from the Broad Institute GATK documenation and involves re-adding the unaligned/unmapped reads into an aligned BAM,  using Picard's MergeBamAlignment.
The BAM files are produced in the end contain both the mapped and the unmapped reads. These files are then sorted with SortSam [2]- so that the sort order in the header becomes:
@HD VN:1.6 SO:coordinate
How does MarkDuplicates handle the unmapped reads of a BAM file containing both unmapped and mapped?
Note MarkDuplicates seems to normally take the BAM's ordering into account, namely, it accepts arguments such as --ASSUME_SORT_ORDER X.  However it's not specified whether reads without a position are ignored, or have to be compared with all other possible reads.

Disclaimer: I initially posted this question on the GATK forum [3], but I'm reaching out to hopefully a broader audience.
Citations:

1 How to Map and clean up short read sequence data efficiently
2 SortSam (Picard)
3 Original posting on the GATK forum

From the Picard documentation:

DUPLICATION METRICS:
Metrics that are calculated during the process of marking duplicates within a stream of SAMRecords.

UNMAPPED_READS    The total number of unmapped reads examined. (Primary, non-supplemental)

It won't alter the flags on these reads, but it will count them in the summary report it generates.
You should be able to test this yourself with a small set of mapped + unmapped reads


swap genotypes for specific fields of a vcf file Is there a tool already established to "swap" the ref/alt information in the GT, AD and PL fields of a vcf file?

For the GT field, I need to transform: 0/0 -> 1/1 ;  1/1 -> 0/0 ;  0/1 -> 1/0 ; and finally 1/0 -> 0/1

For the AD field, I need to swap the columns: x,y -> y,x

For the PL field, I need to swap the columns: x,y,z -> z,y,x

Sample input:

Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  0/0:17,0:17:48:0,48,720 0/1:14,8:22:99:167,0,347 .....

Sample output:

    Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:720,48,0 1/0:8,14:22:99:347,0,167 .....

note: the same operations will need to be carried out for all samples if the VCF file has multiple samples (the ... stand for other potential samples).

Before I dive into scripting ( basically playing with regex), I would like to know if something to process this is already available. I looked into bcftools and fixref but they don't seem to help for this.

If we can always assume your fields will be GT:AD:DP:GQ:PL, you can use this relatively simple, brute approach awk script:

awk -vOFS="\t" -F"\t" '{
  ## just print the header lines
  if(/^#/){
    print;
  }
  else{
    ## The number of samples will be the extra number of fields
    ## in addition to the default VCF 9
    sampleNum = NF - 9
    ## Fix each sample
    for(i=1; i<=sampleNum; i++){
      ## get the fields
      num=split($(9+i), fields,":");
      ## genotype
      if(fields[1] ~ "0/0"){ fields[1] = "1/1"}
      else if(fields[1] == "1/1"){ fields[1] = "0/0" }
      else if(fields[1] == "0/1"){ fields[1] = "1/0" }
      else if(fields[1] == "1/0"){ fields[1] = "0/1" }
      ## AD
      split(fields[2],values,",")
      fields[2] = values[2]","values[1]
      ## PL
      split(fields[2],values,",")
      ## Now set these as the sample data
      $(9+i) = fields[1];
      for(k=2;k<=num;k++){
        $(9+i) = $(9+i)":"fields[k]
      }
    }
    print;
  }
}' file.vcf

On your example input, this returns:

Contig111   11  Contig111_11_T_G    T   G   4278.22 PASS    AC=8;AF=0.073;AN=110;BaseQRankSum=1.35;DP=1176;ExcessHet=0.1448;FS=1.841;InbreedingCoeff=0.1904;QD=13.41;ReadPosRankSum=0.351;set=Intersection  GT:AD:DP:GQ:PL  1/1:0,17:17:48:0,48,720 1/0:8,14:22:99:167,0,347    .....


How to convert mzXML to mzML format (mass spectrometry files)? How to convert mzXML to mzML format (mass spectrometry files) preferably under Linux and/or Python?

ProteoWizard offers this functionality:

msconvert your_mzXML_file --mzML

FileConverter of OpenMS (I believe using the integrated ProteoWizard) also does this.


How to seperate the genus into different group based on th ANI calculation? I've calculated the ANI of the streptomyces genus, but I'm confused about how to separate them into a different group, could anybody kindly tell me how to do that? Thanks!!!

ANI = average nucleotide identity.

This will give you a 2x2 matrix. You reconstruct a phylogeny using clustering, the most important being neighbor-joining, UPGMA is very occassionally used.

The statistic is clustering algorithm and its done by pairwise reconstruction of sister groups. Once a sister-group is defined (usually on nearest-neighor both in neighbor-joining and UPGMA) the genetic divergence from the hypothetical node is established either using the mean (UPGMA) or additive criteria (neighbor-joining). Okay enough theory ...

You get Mega X and find out the import format for a 2x2 distance matric. Mega X however will calculate this for you. You then go to the build tree menu heading and click "neighbor-joining" (or "distance"). The tree will very quickly automatically appear.


Truncating branch length values of Phylogenetic tree with biopython I have been using biopython 1.72 to display my phylogenetic tree files.

Using the function 'Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)' to display branch lengths as well on tree, the tree displays the branch lengths as including all numbers after decimal. While I want the numbers to be truncated to two decimal places.
Following is a python method that truncates upto "the number required" decimal places.

def truncate(n, decimals=0):
    multiplier = 10 ** decimals
    return int(n * multiplier) / multiplier
truncate (12.567, 2)  #Using function
12.56                 #Output
Can I use such method with biopython for branch lengths?
'Phylo.draw' function opens the window which is very narrow and if the tree file is larger, all the branches are overlapped and it gets difficult to read what is written. Does anyone know any better display method?

Following is the example code, I have been using:

from Bio import Phylo
from Bio.Phylo.TreeConstruction import *
from Bio import AlignIO

aln = AlignIO.read(open('example.phy'), 'phylip')

calculator = DistanceCalculator()
dm = calculator.get_distance(aln)

constructor = DistanceTreeConstructor()
njtree = constructor.nj(dm)
starting_tree = njtree
scorer = ParsimonyScorer()
searcher = NNITreeSearcher(scorer)
constructor = ParsimonyTreeConstructor(searcher, starting_tree)
pars_tree = constructor.build_tree(aln)
Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)

I have an msa file as input and it gives a tree, which looks like:

                                                         _______ Human
                                                 _______|
                                                |       |_______ Chimpanzee
                                         _______|
                                        |       |        _______ Dog
                                 _______|       |_______|
                                |       |               |_______ Cow
                         _______|       |
                        |       |       |_______________________ Elephant
                 _______|       |
                |       |       |_______________________________ Mouse
                |       |
         _______|       |_______________________________________ Platypus
        |       |
        |       |        _______________________________________ Anole_lizard
        |       |_______|
  ______|               |        _______________________________ Chicken
 |      |               |_______|
 |      |                       |_______________________________ Zebra_finch
 |      |
 |      |_______________________________________________________ Xenopus
_|
 |       _______________________________________________________ Zebrafish
 |______|
 |      |_______________________________________________________ Fugu
 |
 |______________________________________________________________ Amphioxus

Apologies, try the following,

allclades = list(pars_tree.find_clades(order='level'))
for i, blength in enumerate(allclades):
        printf(blength=%.2f, blength)

"allclades" is a dictionary rather than a simple array.

If this fails I would examine the allclades data structure either using pprint (directly) or using printStruct via the following code,

>>> d = [{'a1':1, 'a2':2, 'a3':3}, [1,2,3], [{'b1':1, 'b2':2}, {'c1':1}], 'd1', 'd2', 'd3']
>>> printStruct(d)

def printStruct(struc, indent=0):
  if isinstance(struc, dict):
    print '  '*indent+'{'
    for key,val in struc.iteritems():
      if isinstance(val, (dict, list, tuple)):
        print '  '*(indent+1) + str(key) + '=> '
        printStruct(val, indent+2)
      else:
        print '  '*(indent+1) + str(key) + '=> ' + str(val)
    print '  '*indent+'}'
  elif isinstance(struc, list):
    print '  '*indent + '['
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ']'
  elif isinstance(struc, tuple):
    print '  '*indent + '('
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ')'
  else: print '  '*indent + str(struc)

Once you know the structure you unravel it accordingly to get to the branch lengths and then issue a printf command to round them to 2 dp.

Goodluck


What do negative dN or dS mean in cal_dn_ds in Biopython? I'm noobie in bioinformatics. I use cal_dn_ds from Bio.codonalign.codonseq with method 'NG86'. Sometimes I get dS=-1. What does it mean? How can number of substitutions be negative? Does it mean that something is wrong with my data?

Upd.
1) Chris_Rands pointed to sources, I read them already, but the question remain. What is that ps variable, and why if it's greater than 3/4, dS = -1?

2) Some more details.
cal_dn_ds requires sequences from CodonAlignment object, which is created by function Bio.codonalign.build that requires (quote from docs) protein alignment and corresponding nucleotides. I have nucleotide sequences from NCBI at the beginning. So I trim non-coding regions from nucleotide sequences, translate all trimmed sequences, write them to a file, perform multiple alignment of the protein sequences, and then pass the list of aligned protein sequences and the list of nucleotide sequences to the build function.
I analize protein DKK1, and example sequences that get negatives are Empidonax traillii XM_027898752.1 and Physeter catodon XM_007110600.1. But I just want to understand general idea behind dS = -1, what can it mean it different cases.

3) Michael G. gave an interesting suggestion that it means zero, but it's strange because the function returns dN and dS separately, so Biopython might as well return zero, so the user would handle those cases himself. And for some sequences Biopython actually returns normal positive dS and dN = 0 at the same time. If it's an error, what are the possible causes? And why in those error cases Biopython can calculate dN, but dS is -1?

dS = -1
.. there are two suggestions

I suggest its an error, possibly based on nucleotide calls within the S, check for anomylous nucleotides, particularly over a non-synonymous codon
It might mean their are no synonymous mutations. It is not the best idea to assign this a value of 0 in a ratio because it will cause the ratio to fall over when in fact their could be strong selection pressure, ie non-synon mutations. Having an incrementally small value is problematic because dN if positive could anything up to infinity, so -1 might flag this problem.

Okay I now know why... You are correct to be concerned with the result because its serious.

My sincere advice to you is do not use any dN/dS ratio where -1 is occurring under any circumstances. On this occasion I'm not going to explain why for two reasons:

First and foremost, as a phylogeneticist you should know why, it relates to the fundamentals of the discipline.
(Hint) The answer is in my past posts on this forum.

To conclude, I'm pretty nuanced in statistical dogma, you can check my record, IMO most dogma is tribalism because every calculation has weaknesses. HERE HOWEVER, there is a real problem and again I urge you to ditch any data set where -1 is occurring and everything associated with it. That is a fairly dire warning should give you a clue about whats wrong. If you are not prepared to do this, then you must report the result and the situations where it is arising (don't hide it).


Convert R RNA-seq data object to a Python object I have done some work in R and would like to try a Python tool.

What is a good way to import the data (and its annotations etc) as a Python object?

I am particularly interested in converting a Seurat object into an AnnData object. (Either directly or as a Python object which can be converted into an AnnData.)

A simple solution for converting Seurat objects into AnnData, as described in this vignette:

library(reticulate)
seuratobject_ad <- Convert(from=seuratobject, to="anndata", filename="seuratobject.h5ad")

Alternatively, one can use Loom, "a file format designed to store and work with single-cell RNA-seq data".

In R, install LoomR:

devtools::install_github(repo = "mojaveazure/loomR")

Convert from Seurat object to loom:

pfile <- Convert(from = pbmc_small, to = "loom", filename = "pbmc_small.loom")
pfile$close() # close loom objects when finished using them.

Then import loom object in Python using loompy, or directly as AnnData:

scanpy.api.read_loom

Alternatively, see feather.

Or export as some text format (csv, json) then import into Python.


Error rate setting in Canu error correction I want to use Canu to correct my nanopore long read (version: MinION R9.5), but I am not quite sure how to set the correctErrorRate. Should I follow the Canu manual (Nanopore R7 2D and Nanopore R9 1D Increase the maximum allowed difference in overlaps from the default of 14.4% to 22.5% with correctedErrorRate=0.225), or you guys have a better option?

It's generally a good idea to trust the "official" suggestions. You can also adjust the error rate based on coverage according to parameter reference:

  For low coverage datasets (less than 30X), we recommend increasing
  correctedErrorRate slightly, by 1% or so.

  For high-coverage datasets (more than 60X), we recommend decreasing
  correctedErrorRate slighly, by 1% or so.

  Raising the correctedErrorRate will increase run time. Likewise,
  decreasing correctedErrorRate will decrease run time, at the risk of
  missing overlaps and fracturing the assembly.

That being said, the sequencing quality can vary from between libraries or flow cells. The error rate is not a constant.

According to the developer (full thread):

  if you want to assemble 1D data or have a bad run, it is possible you
  would need to increase the error rate. You could run with a high error
  rate (0.1) and look at the distribution of overlap error rates in the
  unitig step to look for a peak in the distribution. If you have a near
  neighbor you could also map the corrected reads to it and estimate the
  residual error that way.


Old versions of the reference genome and dbSNP For benchmark purposes I’m trying to find an old version of the human reference genome (pre-GRCh37) and a matching version of dbSNP.

Unfortunately it appears as though older versions of dbSNP aren’t archived, or at least these archives are not publicly accessible.

In particular, the NCBI FTP archive is fairly sparse. The only pre-GRCh37 version it contains is ostensibly dbSNP build 124 from January 2005, for hg17 (NCBI build 35). Unfortunately the corresponding FTP folder is actually empty. dbSNP build 138 has been “removed to recover storage”, whatever that means, and I’m unclear which version of the human reference genome it refers to (its official page does not contain this rather crucial bit of information).

Is there another authoritative source of this data?

The UCSC Genome Browser downloads page has old human genome assembly versions way back to the first human genome ever published (by UCSC in 2000).
http://hgdownload.soe.ucsc.edu/downloads.html#human

We also have a variety of archived dbSNP data up to dbSNP141 on hg38, dbSNP 131 on hg37, dbSNP126 on hg18:
http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/snp141.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/snp138.txt.gz
http://hgdownload.soe.ucsc.edu/goldenPath/hg18/database/snp126.txt.gz

They're in a custom table format, which can be turned into VCF manually.

You can also get all of this data from Table Browser and add filters, etc.:
http://genome.ucsc.edu/cgi-bin/hgTables


How to translate amino acid sequences to Nucleotide sequences I want to convert a list of fasta ( protein sequences) in a .text file into corresponding nucleotide sequences. A Google search gives me result of DNA to protein conversion but not vice versa. Also, I came across How do I find the nucleotide sequence of a protein using Biopython?, but this is what I am not looking for. Is there any possible way to do it using python.Moreover, I would like to solve it using python programming. I am sure there must be some way to do it rather than writing a code from scratch. Thanks!

Although there is not a unique nucleotide sequence that translates to a given protein, one can list all the possible DNA sequences that do translate to that protein.

An online tool that does just that is Backtranambig, from EMBOSS. It produces a DNA sequence representing all the nucleotide sequences matching the input protein, using IUPAC ambiguity codes.


Changing this code in a way to work for my files I have a bash script which extracts some information from a .vcf file but one .vcf file at each time. How I can change this script in a way to work on a bunch of .vcf files and the output is a .txt file named by the name of original .vcf

This is my script

help please

  EDITED

 #!/usr/bash

#outfilename
outname=$".parsed.txt"

#Header for output file
echo -e "$outname\t$chrom\t$Pos\t$Ref\t$Alt\t$ReadCount\t$VariantAlleleCount\t$ReferenceAlleleCount\t$ReadCountControl\t$VariantAlleleCountControl\t$ReferenceAlleleCountControl\t$VAF" >> $outname ;
while read -r line ;
do

#Basic information
chrom=$(echo "$line" | sed 's/ /\t/g' | cut -f 1) #&& echo $chrom;
Pos=$(echo "$line" | sed 's/ /\t/g' |  cut -f 2) #&& echo $Pos;
Ref=$(echo "$line" | sed 's/ /\t/g' | cut -f 4)
Alt=$(echo "$line" | sed 's/ /\t/g' | cut -f 5)

#Tumor sample read, variant and reference information
ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
ReferenceAlleleCount=$(echo $ line | awk -v rc="$ReadCount" -v vac="$VariantAlleleCount" '{print rc-vac}')

#Control or Normal read, variant, reference information
ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

ReferenceAlleleCountControl=$(echo "$line" | awk -v rcc="$ReadCountControl" -v vacc="$VariantAlleleCountControl" '{print rcc-vacc}')

VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
#Print output
echo -e "$outname\t""$chrom"'\t'"$Pos"'\t'"$Ref"'\t'"$Alt"'\t'"$ReadCount"'\t'"$VariantAlleleCount"'\t'"$ReferenceAlleleCount"'\t'"$ReadCountControl"'\t'"$VariantAlleleCountControl"'\t'"$ReferenceAlleleCountControl"'\t'"$VAF" >> $outname ;

#Remove info tags from VCF
done; < <( egrep -v '#' "$1")

And the error I obtained

[fi1d18@cyan01 new]$ bash vcf_parasing.sh
vcf_parasing.sh: line 7: syntax error near unexpected token `echo'
vcf_parasing.sh: line 7: `echo -e $outname'\t' "Chrom"'\t'"Position"'\t'"Ref"'\t'"Alt"'\t'"TumorReadCount"'\t'"TumorVariantAlleleCount"'\t'"TumorReferenceAlleleCount"'\t'"NormalReadCount"'\t'"NormalVariantAlleleCount"'\t'"NormalReferenceAlleleCount"'\t'"VAF" > $outname'
[fi1d18@cyan01 new]$

My script

https://www.dropbox.com/s/di49mij2a43i43s/vcf_parasing.sh?dl=0

Two samples of .vcf

https://www.dropbox.com/s/ryk7q0aybeo5cy2/LP6008202-DNA_B03_vs_LP6008201-DNA_B03.snp.pass.vcf?dl=0

https://www.dropbox.com/s/u44s6d54r5d57bu/LP6008336-DNA_H02_vs_LP6008333-DNA_H02.snp.pass.vcf?dl=0

Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible.

Here's a Perl script that should do exactly what your bash script does but will be several times faster and can deal with multiple input files:

#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile) or
      die("Couldn't open '$outFile' for writing: $!");
  open(my $inFileHandle, '<', $inFile) or
      die("Couldn't open '$inFile' for writing: $!");

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    my @infoFields = split(/;/, $info);
    my $readCount = $infoFields[12];
    $readCount =~ s/ReadCount=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCount = $infoFields[25];
    $variantAlleleCount =~ s/VariantAlleleCount=//;
    #/#Useless comment to fix SE syntax highlighting
    my $referenceAlleleCount= $readCount - $variantAlleleCount;

    my $readCountControl = $infoFields[13];
    $readCountControl =~ s/ReadCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCountControl = $infoFields[26];
    $variantAlleleCountControl =~ s/VariantAlleleCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $referenceAlleleCountControl = $readCountControl - $variantAlleleCountControl;
    my $VAF = $infoFields[27];
    $VAF =~ s/VariantAlleleFrequency=//;
    #/#Useless comment to fix SE syntax highlighting

    my $outString = join("\t", $outFileName, $chrom, $pos, $ref, $alt, $readCount,
                         $variantAlleleCount, $referenceAlleleCount, $readCountControl,
                         $variantAlleleCountControl, $referenceAlleleCountControl, $VAF);
    print $outFileHandle "$outString\n";
  }
  close($outFileHandle);
  close($inFileHandle);
}

Save that file as parser.pl and then run it on all your vcf files with:

perl parser.pl /path/to/vcfs/*vcf

That should create a .parsed.txt file for each input file in the target directory.

If you just have to use the shell for this, you can try with this improved version of your script:

#!/usr/bin/env bash

for file in "$@"; do
  outFile="${file}.parsed.txt"
  outFileName=$(basename "$outFile")
  (
    printf '%s\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t' "$outFileName"
    printf 'TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount\t'
    printf 'NormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n'
  ) > "$outFile"

  echo "Reading $file"
  grep -v '^#' "$file" |
    while read -r line; do
      read chrom Pos Ref Alt rest < <(awk -F"\t" '{print $1,$2,$4,$5}' <<<"$line")
      ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
      VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
      ReferenceAlleleCount=$(($ReadCount - $VariantAlleleCount))

      ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
      VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

      ReferenceAlleleCountControl=$(($ReadCountControl - $VariantAlleleCountControl))

      VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
      printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$outFileName" "$chrom" "$Pos" "$Ref" "$Alt" \
             "$ReadCount" "$VariantAlleleCount" "$ReferenceAlleleCount" "$ReadCountControl" "$VariantAlleleCountControl" \
             "$ReferenceAlleleCountControl" "$VAF" >> "$outFile"
    done
done

Just to illustrate what I mean about the shell being slow, I ran both scripts above on the example files you provided:

$ time foo.pl *vcf 2>/dev/null

real    0m0.051s
user    0m0.047s
sys     0m0.003s

$ time bar.sh *vcf 2>/dev/null

real    2m28.842s
user    3m33.648s
sys     1m12.456s

As you can see above, the perl script took less than a second, while the bash script took almost two and a half minutes!

Finally, if you find the perl script too long, here's the same thing in a more condensed style:

#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile);
  open(my $inFileHandle, '<', $inFile);

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    $info=~s/[^;]+=//g;
    my @infoFields = split(/;/, $info);
    print $outFileHandle (join "\t", $outFileName, $chrom, $pos, $ref, $alt,
                               $infoFields[12],$infoFields[25],
                               $infoFields[12]-$infoFields[25],
                               $infoFields[13],$infoFields[26],
                               $infoFields[13]-$infoFields[26],
                               $infoFields[27]) . "\n";
  }
}


protein secondary structure matching/comparison scoring I want to compare two secondary structures of aligned proteins. I do not know exactly how to do it well.

Example:

EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC

EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHEEEECCCCCCCCCCEEEEEEEEEEEEEE

Example of matching:

EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
|||        |||||| |||                    |||||||||     ||||||      |   |||    ||
EEEHHHHHHHTTTTTTTTCTTHHHHHHHHHHHHHHHHHHHHHHHHHHHHHH----HHEEEECCCCCCCCCCEEEEEEEEEE

Hypothetically: 50% aligned - > score 0.5 - >  thats scoring is wrong:

EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHHHHHHHHHCHCCTHHEEEETTTEEECTTTEEECTTTEEC
||||||||||||||||||||||||||||||||||||||||||
EEEEEEEEEEETTTTTTHCTTTEEECTTTEEECTTECTCCCHEEEEEEEEEEEEEEEHHHHHHHHHHHHHHHHHHHHHHHH

This case will have similar score...

This is a puzzle for me, I don't know how to approach it technically.

Hey i solve this simply by use of SOV alghoritm. https://www.researchgate.net/publication/319526145_A_modified_definition_of_Sov_a_segment-based_measure_for_protein_secondary_structure_prediction_assessment

cheers


I like to use the join command for this. For your example above, you can simply use the following:

$ join file2 file1
ENST00000000412 0.779142 M6PR
ENST00000001008 0.738143 FKBP4
ENST00000002501 0.715315 DBNDD1
ENST00000002596 0.713664 HS3ST1

This assumes your files are sorted. If they aren't sorted, you can use input redirection if your shell supports it (bash and other modern shells do):

join <(sort -k 1b,1 file2) <(sort -k 1b,1 file1)

If your shell doesn't support <(), sort the files first:

sort file1 > file1.sorted; sort file2 > file2.sorted
join <(sort -k 1b,1 file2.sorted) <(sort -k 1b,1 file1.sorted)


How to run enrichment analysis of protein functional annotation? I have a lot of protein clusters. I want to perform an enrichment analysis of their functional annotations, against reference datasets or list of genes I select.

More precisely: a method yields cluster of proteins. I want to decide, for each cluster (which is a set of proteins identifiers), if it holds meaning regarding the proteins it contains, by comparing the annotations found in the cluster and the annotations found in common or specific datasets.

Initially, I used DAVID, which compute the GO annotations from the protein list, then perform the enrichment analysis against common datasets.

However, DAVID suffer of the following drawbacks:

since I didn't find any other way to use it, I'm limited by the web interface and its limitations (number of genes, url size, number of request/day).
automation over the web interface is a hacky script.

The team behind seems to offer a way to run DAVID in-house (allowing, hopefully, to relax the limitations), but I can't find any way to download the database from their website.

What are the alternative ways to get enrichment analysis over proteins in a reliable and automated way ?

I would recommend R instead of sticking to websites.

There are many tools for enrichment analysis. I like to use GOseq, which is initially made for RNAseq data, but can also be used for proteins (I have used it for proteomics). You can use the length of the protein instead of transcript length (to correct for length bias), or you can exclude the length information by just using the hypergeometric method.

See manual for code.

Subsequently you can use gogadget, a tool I made for visualization of GOseq results with a heatmap or cytoscape network.

Additionally, you might wanna take a look at clusterProfiler. With this package you can do enrichment analysis for many clusters in parallel.


How I can test my hypothesises computationally I have single cell RNA-seq data on about 2000 cells in 9 time point. I have clustered my cells in each time point by Seurat. I am seeing in some time points I have 3 clusters while in another time points I have 2 clusters of cells (please look at the plot).
One hypothesis would be that, a new cell type is being emerged at last time point. But maybe if we had another further time point this third cluster being disappeared because that has been arose just because of undifferentiated cells from the later time point (second hypothesis).
Another hypothesis would be this new third cluster is a new cells type that only temporarily expresses a set of genes but if we had cells for another further time point, we would not see the expression of these genes anymore and we had again 2 clusters. I am really struggling with this how I could test my hypothesises.
The final hypothesis would be that the third cluster is really a new cell type born from old cell types during the differentiation.
Like UAD algorithm addresses my question. I just tried FateID algorithm on my data but got stopped by
> pr <- prcurve(y,fb,dr,k=2,m="tsne",trthr=0.4,start=2)
Error in svd(xstar) : a dimension is zero

Any suggestion would be highly appreciated.

This is what packages like FateID and Monocle are for, namely taking single-cell RNAseq data and inferring differentiation trajectories from it. Don't try to reinvent the wheel on this, there are a number of packages out there to do this sort of thing and they're going to have pretty complicated methods that you're really not going to want to reinvent unless you absolutely have to.


Interpretion of my coronavirus 2019-nCov, Wuhan, China BLAST tree? This is the BLAST tree of the latest coronavirus out of China (from Wuhan Institute of Virology, China). It seems strange that there is so much divergence from all the other coronaviruses. Is this expected of new diseases?

The  evolutionary related group (clade) of betacoronaviruses you have identified share an amino acid homology of 85% and include SARS. I know this from the underlying tree published on BioRxiv of a broader group of betacoronaviruses, i.e. your data is a defined subset of the betacoronaviruses which all share a unique, single common ancester.
Lets call this group the SARS clade.
You have performed a nucleotide blast and asked NCBI to produce a NJ tree using 2019-nCov as the reference. I can tell this by the colour coding and the scale bar at the bottom right hand corner shows the genetic distance involved is reasonably larger than 15% divergence by amino acid data. The scale bar represents the number of mutations per nucleotide.
In summary, your tree is essentially a subset of the broader genetic diversity of the one published, but there is a rooting issue.
In your tree, the majority of sequences are from the 2002 SARS epidemic and the virtually zero genetic distance between them is simply because it is a rapidly transmitted outbreak. I didn't realize that SARS has two independent origins both initially from bats. This is quite scary.

2019-nCov is an outgroup within the 'SARS clade' hence it appears at the other side of the tree, i.e. they share a more distant common ancestor
However 2019-nCov is not the most distant common ancester within this subset of the betacoronaviruses, i.e. SARS clade, this belongs to the two bat viruses ZC45 and ZX21. The program has made a likely rooting error (below).
Again the blast omitted the majority of the betacoronaviruses, for example MERS

Rooting issue the reason I suggest there is a rooting error is because the BioRxiv tree using a broad sample of the betacoronaviruses placed the bat strains ZC45 and ZX21 as outgroups to the SARS clade and 2019-nCOV was immediately within that, so in this definition 2019-nCOV is an 'ingroup' within the SARS clade, whereas your tree it is an 'outgroup'. It not a huge issue, but the location of the "root" is determined by the common ancester above it (sequences with <85% homology) and in your tree those have been omitted.
Generally, I like the analysis otherwise and gives an insight into SARS that I hadn't previosly been aware of.


Raw vs Filtered in the output of cellranger count After running cellranger count I got two relevant for further analysis folders: filtered_gene_bc_matrices and raw_gene_bc_matrices. What is the difference between them? What does cellranger filter out, and what should we preferably use for further analysis?

I'm unsure whether this is the answer you are looking for, but when looking into 10X cellranger documentation for the Matrices Output:

  Unfiltered gene-barcode matrices:
  Contains every barcode from fixed list of known-good barcode sequences. This includes background and non-cellular barcodes.

  Filtered gene-barcode matrices:
  Contains only detected cellular barcodes.

Source


Why are some of the gene peptides returned by biomaRt missing an asterisk in the end? Note: this question has also been asked on BioStars

I am getting peptides using biomaRt library in R for the gene 'BRCA1'. 27 different aminoacid sequences are returned, and 12 of those have an asterisk (*) in the end, while 15 do not. What does that mean? I know that asterisk is used to indicate stop codon. So, does that mean the peptides without an asterisk do not have a stop codon? Is that possible at all? Are all of those valid peptides or not?

Below is the R code I used:

library(biomaRt);
mart = useMart('ensembl', dataset='hsapiens_gene_ensembl');
seq = getSequence(id='BRCA1', type='hgnc_symbol', seqType='peptide', mart = mart);

MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKVLLCCPSWSTVVRS* is an example sequence with an asterisk in the end. MDLSALRVEEVQNVINAMQKILECPICLELIKEPVSTKCDHIFCKSLQESTRFSQLVEELLKIICAFQLDTGLEYANSYNFAKKENNSPEHLKDEVSI is an example without an asterisk.

If it's 3' incomplete that means the evidence used to create it was a fragment. Here's the evidence used to construct BRCA1-214 ENST00000477152.5, a 3' incomplete.

You can see that there's a full length cDNA from EMBL, AK307553.1, which was used to create this model.

The sequence was mapped against the genome sequence to create the transcript model. When searching for the translation, it was found that the open reading frame was started, but there was no stop codon. This suggests that the cDNA AK307553.1 is actually a fragment, that the mRNA was broken or cleaved before it was reverse-transcribed sequenced. We display it in Ensembl in the hope that this will lead someone to identify the full length transcript that it represents.

Many people choose not to work with these incomplete transcripts. You can do this by filtering by Gencode Basic in biomaRt.


Odds ratio and enrichment of SNPs in gene regions? I did a QTL analysis with a panel of 7M SNPs, and want to analyze the enrichment of the significant qtl-SNPs in different genic regions (promoters, gene bodies, TFBS, etc.).
A straightforward way to do it would be with an odds ratio (OR), using as 'background' the full 7M SNPs, but I prefer to compare with a background with a similar MAF distribution to that of the qtl-SNPs. For that, I'm thinking in two alternatives:

Take a single random sample of all SNPs with size 10n (n is the number of qtl-SNPs), keeping the same MAF distribution of the qtl-SNPs results, as well as the same proportion of qtl-SNPs in the sample compared to the full 7M panel. Then, test enrichments using ONLY the random sample (with OR and Fisher's exact test).
Do a bootstrapping, taking multiple random samples (excluding qtl-SNPs), each of size n, and keeping as well the same MAF distribution. For each iteration, get the OR between the full qtl-SNPs and the random sample. Then, the mean of the OR will be the point estimate and the quantiles the confidence intervals. For this approach, the problem would be to compute a p-value.

Which approach seems to fit better? If the second one, how to quantify the significance? Can you suggest a better approach?
Thanks.

The first approach only address the question of how likely are you to end up with the observed over-representation given the MAF distribution.
My suggestion is to use the second approach, but I am not sure if you would call it bootstrap. Bootstrap in general means sampling with replacement to estimate the uncertainty of a parameter, in this case, OR. So even if you resample your full dataset, you get an estimate of how variable your enrichment is, not how likely you are to observe it by chance. so just name it permutation or sampling.
This is what I think might work

Calculate your true OR, qtl snps versus all others for say TFBS
Divide your snps according to MAF into 10 bins / stratas
Within each bin, we permute the label whether the snp is a qtl
Calculate OR of selection / others for TFBS
Repeat 1-3 N times to collect the permutated OR
p value would be number of (permutated OR > observed OR+1)/(N+1)

Below is an example implementation in R, first we have 1e5 signficant snps and they are 5 times as likely to be in TFBS:
snp_is_qtl = rep(c(1,0),c(0.1e6,6.9e6))
snp_is_TFBS = c(sample(0:1,0.1e6,replace=TRUE,prob=c(0.95,0.05)),
sample(0:1,6.9e6,replace=TRUE,prob=c(0.99,0.01)))

We have the MAF for all snps, with those in qtl being different from others and put them in bins.
snp_MAF = c(runif(0.1e6,0.4,0.8),runif(6.9e6,0.1,0.9))

We sort the snps according to MAF to make it easier:
o = order(snp_MAF)
snp_is_qtl = snp_is_qtl[o]
snp_is_TFBS = snp_is_TFBS[o]
snp_MAF = snp_MAF[o]

And we create the bin:
bin_MAF = cut(snp_MAF,10)

Lastly a function to calculate OR:
calOR = function(snps,TBFS){
tab = table(snps,TBFS)
log(tab[1,1])+log(tab[2,2]) -log(tab[1,2]) -log(tab[2,1])
}

So our observed OR (in log) is:
obs_OR = calOR(snp_is_qtl,snp_is_TFBS)
obs_OR
[1] 1.656859

Now we collect the OR under permutation:
N = 100
permutated_OR = sapply(1:100,function(i){
print(i)
perm_qtl = unlist(tapply(snp_is_qtl,bin_MAF,sample))
calOR(perm_qtl,snp_is_TFBS)
})

hist(permutated_OR)

Your observed OR way exceeds the permuted, even after controlling for MAF. Note the computation is a bit heavy because of the large number of snps. So it makes sense to calculate all the ORs (i.e TFBS, promoters..) with each permutation


What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line I have two DNA strings:

  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC

and

  AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I want a tool that allows me to do something like this on the command line:

$ aligner  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

and receive an ASCII visualization of a pairwise alignment. Something like this would work:

GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---------
---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I prefer a tool that is packaged on bioconda. A tool that also does protein and RNA sequences is even more preferable.

You are looking for the needle program from the EMBOSS suite. Available in bioconda.

http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html

To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, extend_penalty 0.5) you'll need -auto. Thus your query above would be:

$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
#
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47

#---------------------------------------
#---------------------------------------


Should the cell sorting marker genes be excluded during clustering? We sort different populations of blood cells using a number of fluorescent flow cytometry markers and then sequence RNA. We want to see what the transcriptome tells us about the similarity and relation between these cells. In my experience on bulk RNA-seq data, there is a very good agreement between flow cytometry and mRNA expression for the markers. It's good to remember we sort using a few markers only, while there are hundreds if not thousands of cell surface proteins.

Should we exclude genes of these sorting (CD) marker proteins when we perform PCA or other types of clustering?

The argument is that after sorting, these genes may dominate clustering results, even if the rest of the transcriptome would tell otherwise, and thus falsely confirm similarity relations inferred from sorting.

Should we at least check the influence of these genes on clustering results?

I do not think there is a simple "yes" or "no" answer here.

A good starting point would be, as you suggest, use all the genes and assess the results in the light of the marker genes and expected results. This could both serve as as good quality control as well as give you overview of all the processes happening in the cells.

Depending on the marker genes effect and biological question you may then want to remove the marker genes with potential ordering effect, or even other genes, e.g. by GO terms or pathways.

You most likely want to account for the cell cycle phases as well. And check for other technical factors. I recommend Bioconductor pipeline to get inspired https://www.bioconductor.org/help/workflows/simpleSingleCell/ when it comes to scRNA-seq analyses.


Optitype for Singularity I am trying to use optitype on an HPC without docker support or sudo privileges. A search brought me https://github.com/researchapps/optitype and https://github.com/nf-core/hlatyping. I tried bootstrapping the singularity image as follows:
singularity build optitype.simg docker://fred2/optitype

However, this does not seem to run as expected. Could someone with experience bootstrapping docker files onto singularity suggest how I could change this? How do the commands change when doing a bootstrap? I would appreciate your help!

Bit of an old post sorry, but I think you might be getting a little confused here. The nf-core pipelines come with their own software already, so you don't need to install or bootstrap optitype yourself. If it's used in the pipeline, it comes in the pipeline container(s) (or conda environment or whatever packaging system you use). So you should never have to touch docker://fred2/optitype.
When you run the pipeline with the -profile singularity argument, it will automatically pull the image from Docker Hub and convert to Singularity for you (some of the newer pipelines pull pre-build Singularity images directly). You can also use the nf-core download hlatyping --singularity command from the helper tools to pull the image ahead of time.

Docker images for this pipeline: https://hub.docker.com/r/nfcore/hlatyping
Documentation for this pipeline: https://nf-co.re/hlatyping/
nf-core/tools download docs: https://nf-co.re/tools/#downloading-pipelines-for-offline-use

I hope this helps! Come and join the nf-core Slack organisation if you need more help. We have a channel for each pipeline, including #hlatyping.


STRINGdb: reverse map STRING identifiers to gene names? Is there a functionality within STRINGdb package or in another environment to reverse map the STRING identifiers to Gene names?
For example, the entire STRING network can be downloaded using the commands:
library(STRINGdb)
string_db <- STRINGdb$new(version="11", species=9606,
                          score_threshold=0, input_directory="" )
full.graph <- string_db$load_all()
netFull
# IGRAPH c07d6e9 UN-- 19354 11759454 --
# + attr: name (v/c), combined_score (e/n)
# + edges from c07d6e9 (vertex names):
#  [1] 9606.ENSP00000000233--9606.ENSP00000272298 9606.ENSP00000000233--9606.ENSP00000253401
#  [3] 9606.ENSP00000000233--9606.ENSP00000401445 9606.ENSP00000000233--9606.ENSP00000418915
#  [5] 9606.ENSP00000000233--9606.ENSP00000327801 9606.ENSP00000000233--9606.ENSP00000466298
#  [7] 9606.ENSP00000000233--9606.ENSP00000232564 9606.ENSP00000000233--9606.ENSP00000393379
#  [9] 9606.ENSP00000000233--9606.ENSP00000371253 9606.ENSP00000000233--9606.ENSP00000373713
# [11] 9606.ENSP00000000233--9606.ENSP00000248929 9606.ENSP00000000233--9606.ENSP00000356737
# [13] 9606.ENSP00000000233--9606.ENSP00000338766 9606.ENSP00000000233--9606.ENSP00000386135
# [15] 9606.ENSP00000000233--9606.ENSP00000319192 9606.ENSP00000000233--9606.ENSP00000432568
# + ... omitted several edges

V(netFull)$name will provide all the STRING identifiers within this network netFull. How to get the corresponding gene names?
For example, string_db$map() maps the gene identifiers of the input dataframe to STRING identifiers. It returns the input dataframe with the "STRING_id" as an additional column. Now, I would like to explore in the reverse direction going from STRING identifiers to gene names.

This answer does not rely on STRINGdb package.
STRING provides supplementary identifier mapping files here: https://string-db.org/mapping_files/.
It looks like you are interested in human STRING name one - https://string-db.org/mapping_files/STRING_display_names/human.name_2_string.tsv.gz. You could read it as a separate lookup.


ComBat for batch effects removal on scRNA-seq data Is it possible to use sva's ComBat for batch effects removal on scRNA-seq data?

Is there any difference between RNA-seq and scRNA-seq data which doesn't allow to use ComBat on single-cell data? (I am thinking about extreme sparsity, > 80%). If so, are there alternatives?

Look at this recent paper that uses ComBat on scRNA-seq data for batch effect removal and states that it "successfully does so".

I also suggest that you check out this publication on Distribution Matching Residual-Nets. Authors evaluated their method also on scRNA-seq data and thus it may be something you are looking for. I personally played a bit with their code and found it straightforward to use.