Proteomics data Vs Transcriptomics data? I want to use either of Proteomics or Transcriptomics data for integrating it into my kinetic model. Before proceeding, I want to know what are the advantages of using either of them so that I could make an informed decision on it!

Many studies have shown that the best we can do is integrating both transcriptomics and proteomics data with our kinetic model, but I've some time constraints and have to proceed with only one of those.

My effort and findings: I've found from discussions with researchers that gathering transcriptomics data has an amplification step which increases the chance of finding a particular one whereas gathering proteomics data has no such step but has fragmentation and then rejoining which creates many problems(such as splice variants etc) and thus leads to a loss of data. But a PostDoc told me that even after the loss of that data, I'll get more information from Proteomics data.

I want to know such type of points and want to know if these are valid or not!

This is what I found from doing some research. Comments are welcome at any point!

Capture percentage for data gathering:

Transcriptomics Data: There's an amplification step in Transcriptomics data gathering methods. Hence, it's possible to capture almost the totality of the Transcriptome using those methods(scRNA seq methods, Nanopore tech, Spatial transcriptomics, etc)
Proteomics Data: There's no amplification step in Proteomics data gathering methods. Those methods have fragmentation and a rejoining step which leads to a loss of the data.

Scalability of methods:

Transcriptomics Data: Scalable methods available
Proteomics Data: Less scalable methods for protein studies

Reference availability:

Transcriptomics Data: Fully annotated references available on consortiums such as ensembl biomart, etc
Proteomics Data: Universal and Comprehensive Human Proteome reference is still in question and because of the incomplete and inaccurate references, much of the newly generated data is being rendered useless

Uniformity

Transcriptomics Data: Lots of publications on pipelines, and their benchmarking available
Proteomics Data: Lack of uniformity across labs/research groups and lack of related literature lead to differences in protein fragmentation and solubilization and differences in algorithms to run analyses

Technical Bias

Transcriptomics Data: Many methods available to tackle cell-bias, noise in data, and batch effects, so as to capture maximum biological variability
Proteomics Data: Cell-wide Mass spectrometry has a bias towards the identification of peptides with higher concentration or contamination from other experiments


Two variants associated with same chromosomal position in genotype data for one individual I'm trying to play around with my family's raw genetic data from Ancestry.com so that I can do some genetic analyses, and I noticed in the AncestryDNA raw data files that some regions of the genome have more than one variant associated with them and I would like to know why this is. The end goal is get this into a regularly formatted VCF, but I feel that these double variants will probably mess things up. So for example one might see something similar to this in the tab-delimited file:

rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

(Header line: rsid    chromosome    position    allele1    allele2)

I understand that when rsid numbers were created, sometimes multiple rsid numbers were accidentally assigned to the same region of the genome, but  I do not understand how two separate variants can be assigned to one single region in the genome. I could be wrong but as I understand it, I should only see one line here (either a "T    T" or an "I    I") because genotypes are reported in comparison to the human reference genome and so each position should uniquely refer to a single nucleotide in the said human reference genome.

But in this example case one line says there is a "T    T" genotype at this position, and the other line says there is an insertion(see "How do I read my raw DNA data?" here).  One possible explanation I thought of is if this would just be saying that a T was inserted, but if this is the case, then as I understand it the reference allele present in the individual should still reported, just in conjunction with the allele at the index. (So assuming an extra "A" was inserted after the T, which could be the reference or alternative allele, something like this should be recorded: rs6010164  22  50515270    TA  TA).

Can anyone help me understand why there are two separate variants associated with chromosomal position? Thanks!

It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here.

Without info about the reference allele it's pretty confusing (particularly because they appear to be reporting the opposite strand from the database (G vs. C/T), but I would interpret this to mean that both alleles match the reference for rs267607236 (they have the insertion), but do not match the reference for rs6010164 (G-C pair rather than T-A). This is also the most common genotype at this position, being found in almost 80% of people in most populations.


Get raw genome sequence I'm currently working on testing some matching algorithms for strings. I would like to do some tests on raw genomic data as I expect different results from random strings given the lower entropy. How can I get a whole human genome sequence? (just the raw basis)
I have tried websites like https://www.ncbi.nlm.nih.gov/ it seems very confusing and full of unnecessary information given that I just want to have the raw basis without annotations.

You want FASTA files from here.
If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).
See the README for more info.


Extracting all information about a sample when using xtract from e-utilities I would like to extract all information about each SAMPLE after running the following query (run the query and add a | grep SAMPLE for clarification on what I mean by SAMPLE):
esearch -db sra -query PRJNA514750 | efetch -format xml
I tried the following:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern EXPERIMENT -element SAMPLE
but this returns nothing (PS: SAMPLEs are within an EXPERIMENT tag). I read in the e-utilities guide that -pattern will divide the data into rows, and -element into columns, so I'm presuming that this didn't work because SAMPLE has multiple tags within it. So I then tried:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element random_SAMPLE_tag

where random_SAMPLE_tag is any tag within SAMPLE.
Here's a concrete example:
esearch -db sra -query PRJNA514750 | efetch -format xml | xtract -pattern SAMPLE -element TITLE

This works, but I want to get all the information about each SAMPLE, and I do not know beforehand what the tags within it are (I manually got TITLE in this case), and since I want to get this info for a quite a few studies, I can't manually check this.

You can give pysradb a shot.
study_accession experiment_accession experiment_title                                                   experiment_desc
                                      organism_taxid  organism_name           library_strategy library_source  library_selection library_layout sample_accession sample_title instrument           total_spots total_size  run_accession run_total_spots run_total_bases
SRP178555       SRX5243234             GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq    GSM3560726: 97-gem-Day3-OEA115; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245909                    Illumina HiSeq 2500  19052916     446793895  SRR8435699    19052916         971698716
SRP178555       SRX5243233             GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq    GSM3560725: 97-gem-Day3-OEA087; Caenorhabditis elegans; RNA-Seq  6239            Caenorhabditis elegans  RNA-Seq          TRANSCRIPTOMIC  cDNA              SINGLE         SRS4245908                    Illumina HiSeq 2500  18156224     428457264  SRR8435698    18156224         925967424
.. [truncated]


Regarding canonical analysis of principal coordinates (CAP) for linking microbial data with elements I used the below code to do CAP analysis to link microbial data with elements, but not sure how to do statistical analysis to show which element variation is responsible most for the difference in microbial composition. I can only able to see two axis in figure so how I can get complete information about the variation based on these mineral element and individually like which mineral responsible most for microbial diversity?
    # CAP ordinate
    cap_ord <- ordinate(
        physeq = phyloseq_prop,
        method = "CAP",
        distance = bray_not_na,
        formula = ~ ~ P + Mn + Cu +Zn + Rb + Mg )

        # CAP plot
    cap_plot <- plot_ordination(
      physeq =phyloseq_prop ,
      ordination = cap_ord,
        color = "SampleType",
        axes = c(1,2)
    ) +
        aes(shape = Region) +
        geom_point(aes(colour = SampleType), alpha = 0.4, size = 4) +
        geom_point(colour = "grey90", size = 1.5) +
        scale_color_manual(values = c("#CBD588", "#5F7FC7", "orange","#DA5724", "#508578", "#CD9BCD",
       "#AD6F3B", "#673770","#D14285","#652926", "#C84248",
      "#8569D5")
        )

        # Now add the environmental variables as arrows
    arrowmat <- vegan::scores(cap_ord, display = "bp")

    # Add labels, make a data.frame
    arrowdf <- data.frame(labels = rownames(arrowmat), arrowmat)

    # Define the arrow aesthetic mapping
    arrow_map <- aes(xend = CAP1,
        yend = CAP2,
        x = 0,
        y = 0,
        shape = NULL,
        color = NULL,
        label = labels)

    label_map <- aes(x = 1.3 * CAP1,
        y = 1.3 * CAP2,
        shape = NULL,
        color = NULL,
        label = labels)

    arrowhead = arrow(length = unit(0.02, "npc"))

    pdf("cap2.pdf", height=8, width=8)
    # Make a new graphic
    cap_plot +
      geom_segment(
        mapping = arrow_map,
        size = .5,
        data = arrowdf,
        color = "gray",
        arrow = arrowhead
      ) +
      geom_text(
        mapping = label_map,
        size = 4,
        data = arrowdf,
        show.legend = FALSE
      )

      dev.off();

     cap_ord

     Call: capscale(formula = distance ~  P + Mn + Cu +Zn + Rb + Mg, data = data)

                                  Inertia Proportion Rank
Total          2.069e+01  1.000e+00
Constrained    5.240e+00  2.533e-01   10
Unconstrained  1.545e+01  7.468e-01   97
Imaginary     -1.090e-03 -5.268e-05    1
Inertia is squared Bray distance

Eigenvalues for constrained axes:
  CAP1   CAP2   CAP3   CAP4   CAP5   CAP6   CAP7   CAP8   CAP9  CAP10
2.5733 0.8789 0.5647 0.3797 0.3037 0.1873 0.1243 0.1030 0.0800 0.0453

Eigenvalues for unconstrained axes:
 MDS1  MDS2  MDS3  MDS4  MDS5  MDS6  MDS7  MDS8
4.790 1.343 1.138 0.916 0.784 0.649 0.518 0.501
(Showing 8 of 97 unconstrained eigenvalues)

Many thanks

I don't think that you can say that the 25% can be said to relate to the chemical variable specifically; rather, that 25% is the variance that is explained by all constrained axes. While these axes are drawn in such a way that they maximize variance associated with your independent variables, I do not think that they are directly interpretable as being 100% due to the elements.
The number of constrained axes is a parameter chosen as an input to the CAPscale method:

An extremely important point is to determine how many PCO axes should be retained; i.e., what should be the choice for $m$? One might be concerned that the method could suffer from a poor choice for $m$. If $m$ is too small, then there may be some ecologically important information in the data that will not be included in the canonical analysis. On the other hand, if $m$ is too large relative to $N$, then a misleading canonical plot could result. For example, we cannot set $m = N - 1$, for then the squared canonical correlations would all be equal to 1.0, because we have $N$ points!

So the value of the variance there will vary depending on how many constrained axes you use. If I understand right, you could obtain 100% distance variance explained by choosing as many constrained axes as you have observations.
I think that it is worth looking into how others have solved this problem. For example, one has to first decide whether or not this 25% is an unexpectedly large number, i.e. compute  p-values (I am not crazy about p-values in most contexts but here they could serve as a helpful heuristic), as suggested by the CAPscale paper:

As a further option, one may test the hypothesis of either (1) no significant differences in multivariate location among groups (CDA), or (2) no significant relationship with quantitative environmental or other variables (CCorA). This is done by using the trace statistic (sum of canonical eigenvalues (sum of squared canonical correlations; see the Appendix) and obtaining a P value by permutation (e.g., Anderson 2001b).

This seems to be directly implemented as the anova.cca() function in the vegan package (I believe that they have decided that CAP is a flavor of CCA for that package). I understand the motivation to get % variance as it's intuitive, but in complex statistical analysis of distances it is a little bit less clear what it really means. Some kind of comparison to a null hypothesis is probably necessary to say whether it is meaningful.
For how others have approached similar questions you can see this discussion. For other implementations that may be more user-friendly, you might consider the vegan package, which is a richer toolbox for ordination in ecology (even if it is less focused on metagenomics). Here is the CAPscale function in vegan. The phyloseq authors have correctly decided to build large pieces of phyloseq out of vegan, so I think that you can pass phyloseq back into vegan functions and e.g. use those analyses directly on your cap_ord object. But see that link for more info.
For a little more easy-to-understand background on statistics of distance methods (particularly in microbiomes), you might see here.


How to upload and process a data set to Galaxy instance and retrieve results programmatically? How to upload and process (de novo assembly pipe) a data set to Galaxy instance and retrieve the results programmatically? Are there any packages beyond BioBlend?

You best bet is bioblend, since it's a more convenient wrapper around the Galaxy API, but if you really want to you could use that directly. I should note that large files should probably still be uploaded via ftp or sftp rather than via the http(s) interface.


What does the number mean in an HGVSp annotation? Let's take the example p.Arg452Pro that I got from an annotated VCF file that contains an HGVSp column.

What does 452 mean in this case ?

Theoretically, is it possible to have the same HGVSp on the same gene but on a different transcript level ?

For example:

Chromosome      Position        Ensembl_Gene        Ensembl_transcript      HGVSp
4               4613173         ENSG00052342        ENST00141456            p.Arg452Pro
4               4613176         ENSG00052342        ENST00141458            p.Arg452Pro

Thank you in advance.

It's the position in the protein. p.Arg452Pro means that there is a variant that changes an Arginine residue at position 452 of the protein into a Proline.

And yes, it is absolutely possible to have the same HGVS p. notation on the same gene but on different transcripts. More than possible, it is quite common. Your particular example isn't actually a real variant since that position on chromosome 4 is intronic in both hg19 and hg38, but TP53:p.R175L1 is a real example of what you describe:

As you can see in the image or by visiting the link, this variant on TP53 has the same HGVS p. notation on 6 transcripts. It also has a different p. notation on two more transcripts that are not shown in the image as well as no p. notation for a few other transcripts in which it doesn't affect the protein product.

1Disclaimer: The link goes to https://varsome.com, a free human variant search engine that is developed by the company I work for.


What is PyMOLs search path for scripts? I know I can run my PyMOL scripts if I specify a valid path from the current working directory and also using absolute paths, but can I have a central script repository, so I don't always need to specify the full path or copy my scripts around from project to project?

Environment info:

Kernel: Linux 5.6.8-1-MANJARO x86_64
Python: 3.8.2
Python2: 2.7.18
PyMOL: 2.3.0 Open-Source

You can use pymolrc and create aliases to your scripts.

If your script is located e.g. in

/Users/tillyboy/pymolscripts/myscript1.pml

Just enter this line to your pymolrc:

alias myscript1, run /Users/tillyboy/pymolscripts/script1.pml

Now myscript1 is available for all your Pymol sessions independent of your current working directory.


I have used TER to break the long bonds of a chain in my PDB I’m now not sure what I need to alter in my PDB to get it to work in leap. I know breaking the bonds turns the formerly connected residues into terminal residues.
It keeps saying that 3 of my atoms no longer have a type , how can I fix this ?

Normally you have three options in a PDB:

Add a TER record
Add the  OXT atom, if your chain was cleaved or cleaved spontaneously (i.e. Asp-Pro sequence), the C-terminal would sill have a terminal oxygen. i.e. making it a terminal residue. If you have protons, you'll need to add also 3H .
Assign a different chain

If you are using forcefields as is your case, the second is the only choice. The C  atom and the N  atoms have a type that does not match because they have the wrong number of bonds —hence the type error.


Are phylogenetic tree construction algorithms any different than general clustering algorithms? Are phylogenetic tree construction algorithms any different from clustering algorhithms? I suspect the answer is no.
Of course phylogenetic tree construction uses biological knowledge, e.g special distance metrics, but does it brings anything new to clustering algorithm (e.g hierarchical, neighbor-joining etc).

The goal of a phylogeny is to estimate the "expected" number of mutations between all taxa in the analysis and their hypothetical common ancestors. A cluster-analysis will only identify the "observed" mutations and "expected" and "observed" mutations can be majorly different due to the major artefact of reversion mutation. This is particularly true of nucleotide phylogenies.

The key difference between clustering algorithms based on a "uncorrected" distance matrix and phylogeny is the latter is based on an explicit model to accommodate reversion mutations. The real problem is that there are only 4 bases, so randomly there is 1/4 chance that a mutation at a given position will mutate back to the original, e.g. A->C->A. Observed mutational differences = 0, expected (real) mutational differences = 2. What phylogeny is concerned about is reconstructing that "2". The problem is significant because almost every gene has regions of rapid mutations and regions of low mutations.

The primary way of doing this is via the Jukes-Cantor correction and is foundational in all nucleotide trees expect "p-distances". If nucleotide divergence is lower than 75% then it is possible to estimate the expected number of mutations using the observed via the JC correction. Moreover, when combined with a method of estimating rate variation within a gene, (usually the discrete gamma distribution), JC correction is very effective at recovering the "true tree". This is because fast evolving homologous sites are grouped together - big correction via JC, slow evolving sites are grouped together - small correction via JC. Other approaches to improving the JC correction are via identifying the bias between purine to purine mutations and pyrimidine to pyrimidine and purine to pyrimidine mutations.

The importance of JC was demonstrated by simulation studies of 4 taxa ((a,b), (c,d)), if two strains evolved very quickly when there sister lineages evolve slowly a clustering algorithm will report the fast strains are a sister group, i.e. ((b, c),(a, d)). If the JC method is implemented via maximum likelihood (or Bayesian) then it correctly recovers the true tree ((a,b), (c,d)). The artefact is known as long branch attraction.

Clustering algorithms based on a distance matrix that implements the JC correction, tend to perform poorly for long branch-attraction artefacts. This doesn't mean the correction is useless, just not particularly powerful. The issue is that distance matrix methods don't "stick" to the model and the clustering will introduce a layer of inaccuracy. Normally presentation of a "corrected" distance matrix combined with a clustering algorithm will require bootstrapping (resampling with replacement) to assess whether a given cluster is supported. Parameterised distance matrices, using neighbour-joining clustering in in-conjunction with a bootstrap are considered okay.

The @user172818 mentioned parsimony and this method is considered less reliable because it can't implement a JC correction. IMO it is possible weighted parsimony could make a "come-back", but it would be really complicated to implement a biological weighting method and would require extensive, independent calculations.


Fasterq-dump: --split-spot or -concatenate-reads? After using files that I downloaded from the SRA with fasterq-dump, I realize I am not 100% sure that I have all the data.

I noticed in my downstream analysis that I seem to be missing the .1 and .2 numbers in the code associated with individual reads with the same spot ID. Then I noticed that my initial fastq file doesn't seem to have those either.

So I'm wondering if the -split-spot option that I used is doing what I want it do do, which is 'get all the sequence data', and how does it differ from --concatenate-reads. I am afraid I did not understand this from the program description. Please help?

Fasterq comes from the latest version of sratools. So if you check the manual , it says the equivalence is:

fastq-dump SRRXXXXXX --split-3 --skip-technical
fasterq-dump SRRXXXXXX

In older versions of sratoolkit, if you use fastq-dump without specifying --split-3 for paired-end reads, you get the format mentioned, spotID.1 for forward, spotID.2 for reverse:
fastq-dump SRR6462984.sra
more SRR6462984.fastq

You get:
@SRR6462984.1 1 length=301
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGGATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG
+SRR6462984.1 1 length=301
#/AAAEE/EEEA/EAAEEAEE/EE/EEA/EEEEE/AAE/EE/EAEAE//EEAEE/EEEEE<EEEEE/EE<EE/AEEEE/EAEEEEA/E<EEAEAEEEEEEAEEEE<EEAEEEEEAAAEEE/EEEEE/E//AEEEE6EA66EE<<<EAAA/AAAAAE/E6E/EA/AE/AEEEEEEAAEEEEA6EEEEEEEEE/EEE/AEEAAAEEEAE/E/EEEEEA//E//EEE/EEE/A<//AE/EE//6A///AA/E//E/E/EE//AE//AAEEE////E/A<AAE/EEAE//A//6/A</<6/</E/
@SRR6462984.2 2 length=302
NTCAGCAGGCGCCGTAACTTCAGCAAAACTGTAGCTATAACTTTTCCGAACTTTATACCTAACTTCGCTTTCCTAGCACTGTATAAAATGCCTAAAAGAAATCAGCAGGCGCCGTCACTTCAAAAAATCAGTATGGATCAATCAGACGCGGAGCTTATTAAAACAGAGCTTATTAAAAATTATGTCTTTCTTCTTCTTTCAGCAATTCGCAACGTAAAATTCAAAATTTTGGACGCCGTAATTGCGCAGGAACCAATTCAAAGAGGCGGCGGCCAATTAATTCCAACGTCTCGAAGAGTGGC
+SRR6462984.2 2 length=302

With the latest version, it by defaults splits the paired end sra file into forward and reverse fastq (i.e --split-3):
fasterq-dump SRR6462984.sra
head -2 SRR6462984_*.fastq

And you get
==> SRR6462984_1.fastq <==
@SRR6462984.1 1 length=150
NATCTGCCCGTTCCACATAAACACCGTTAAACATCGATAGAGCGAAATAAAGCTGCTGAGTGATGTGATGATAAAGCTTTTCCGTTCTGACACGATCTTCCATCTCGTCTATCATTGCATCGAGAGCGTCAGAGATCGCAAGCAGTGCGG

==> SRR6462984_2.fastq <==
@SRR6462984.1 1 length=151
ATATTGCCGTCCGTCTTGAAGCTTTGGCAGAGCCAGGCGGAATTTGTTTATCGGATAGTGTTTATGCGCAAATTTAGCATATGATTTCCGACAACTTTTGCTTCATCGGGCCGCAAAGCGTGAAAAACAGTGTCGTCACCGTTGATGTGTG

So if you need to get back the .1 and .2 format, use fastq-dump from the latest installation.
I would suggest sticking to the newest format, with 2 files, because fastq-dump will be deprecated in the future.


Gene expression data in GC-RMA format? I'm having a hard time retrieving the most expressed genes in a cell line. One of the databases I've found has values for a lot of genes in 60 different cell lines and it seems the data is obtained using GC-RMA. The data is available here in the file NCI60_U133A_20070815.raw.csv.

My problem is what does the values actually represent? Would it be correct to assume the genes with high values are the ones most expressed in the corresponding cell line?

I've found a short explanation on why there is no unit definition on the page here.

Here is a snapshot of the file:

-            huh-7    huh-7      HEK293  ...
1007_s_at    31.727    55.551    360.731
1053_at     402.633   441.622    557.782
117_at        9.968    10.178      9.565
...

These values come from a microarray so they represent the number of sequences that where attached to the microarray's probes (1007_s_at, ...).

In general yes, a higher value means more expression. The normalizations (GC-RMA) applied are to correct for several effects. They help for testing the probability of the given probe to be a false positive.


hg38 GTF file with RefSeq annotations I'm not sure what I'm missing, but I'm struggling to find an official hg38 GTF file with RefSeq annotations. I'd like to provide the GTF to Salmon to get gene-level annotations.

Here's Salmon's help info for --geneMap:

  File containing a mapping of
                                          transcripts to genes.  If this file is
                                          provided Salmon will output both
                                          quant.sf and quant.genes.sf files,
                                          where the latter contains aggregated
                                          gene-level abundance estimates.  The
                                          transcript to gene mapping should be
                                          provided as either a GTF file, or a in
                                          a simple tab-delimited format where
                                          each line contains the name of a
                                          transcript and the gene to which it
                                          belongs separated by a tab.  The
                                          extension of the file is used to
                                          determine how the file should be
                                          parsed.  Files ending in '.gtf', '.gff'
                                          or '.gff3' are assumed to be in GTF
                                          format; files with any other extension
                                          are assumed to be in the simple format.
                                          In GTF / GFF format, the
                                          "transcript_id" is assumed to contain
                                          the transcript identifier and the
                                          "gene_id" is assumed to contain the
                                          corresponding gene identifier.

Here's what I've tried:

NCBI seems to only provide GFF files. I know you can convert between the two, but that hasn't worked well for me in the past. Both file formats allow a lot of freedom, which makes conversions sloppy. Salmon, specifically, is looking for the gene_id and transcript_id annotations.
UCSC provides GTF files from RefSeq, but the gene_id annotation is identical to the transcript_id annotation (i.e., it's the NM number). Or maybe there's an option I'm missing.
Ensembl provides exactly the format and all of the information I could hope for, except the transcript_id is the Ensembl ID (ENST#), not RefSeq (NM#).

Can anyone tell me how to get a GTF file with hg38 RefSeq annotations?

You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format

The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.

This is the command I usually give when we get this question on our mailing list:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

General questions: genome@soe.ucsc.edu
Questions involving private data: genome-www@soe.ucsc.edu
Questions involving mirror sites: genome-mirror@ose.ucsc.edu


Hard-clip primers in a bam file I am working on a sample which was sequenced using an amplicon kit and I would like to ensure that no variants are called in the PCR primer sequences. I know that the primer sequences are part of the genome and a variant there is no less valid than a variant anywhere else, but for reasons I cannot control, I need to ensure that no variants can be called in the primer regions.

My first thought was to use a soft-clipping tool like BAMClipper. However, the variant caller I am using (Sentieon's implementation of GATK) preforms a local realignment when calling and if it finds that soft-clipped regions align perfectly against the genome, it will remove the soft clips. This makes perfect sense for cases where the softclips were added by the aligner, but sadly breaks my workflow since I need the softclipped primers to remain clipped.

I would therefore like to hard clip, to completely remove the primer sequences from the bam file before variant calling. I have found one tool that seems to be able to do this, TrimPrimers by Fulcrum Genomics, but this fails on my input bam file with an error I can't really make much of.

Is there any other tool that can take a list of regions (e.g. a bedpe file) and a BAM and hard-clip the regions from the BAM?

Another option would be to adjust the quality of the primer bases to be below the threshold that your variant caller uses.  In my case I have "masked" the primer bases by dropping the quality to 0.  This prevents the variants from being called based that are synthetic but if you have an overlapping amplicon it will call bases from those reads.


Authoritative source on human cytogenetic regions? I am looking for a database that would keep track of  human cytogenetic regions and genomic coordinates per genome assembly. I had expected the Genome Reference Consortium to have it, or Ensembl, or NCBI, but I could not find it.

Is there a place with a translation table between human cytogenetic regions and genomic coordinates?

EDIT

After benn's answer I found that you can also get this information from Ensembl's REST API: https://rest.ensembl.org/info/assembly/homo_sapiens?content-type=application/json&bands=1.

I have transformed the data provided by Ensembl and generated a CSV file with the cytogenetic regions:

https://github.com/ramiromagno/gwasrapidd/blob/master/data-raw/cytogenetic_bands.csv

Do you mean cytoband coordinates? For the human genome you can find them for example at UCSC.

hgdownload.cse.ucsc.edu/goldenPath/hg38/database/cytoBand.txt.gz


Clustering of gene co-expression network by igraph R package I have constructed a gene co-expression network from RNA-seq data. The network has more than 10000 nodes and more than 1 million edges.

The network file as an edge list format of memory around 1gb which was created by calculating Pearson correlation of each gene pairs and gene pairs having correlation greater than 95% we're selected to create the edge list. The file is too large to run in cytoscape due to limitation of computer configuration.

In R, I am able to calculate the topological properties of the network. And to cluster the gene coexperssion network, I tried using different community detection methods and but I am unable to identify which community detection is good for obtaining clusters?.

After clustering how to rank the subclusters and save each subcluster in separate files.

For the first part, do you mean that the file is too large to be run on your computer?

For the second, if I understood correctly, you can use igraph or ggnet2 to color code the nodes based on the conditions you have and see if communities are formed based on condition. This would be a simple way to do this kind of thing.

I personally prefer using ggnet2 as I found it more user friendly. What you can do if you have an adjacency matrix, say adj_mat is to create a color vector based on the conditions for each node in the network, where length(color_vec) = ncol(adj_mat) = nrow(adj_mat) - since it is a square matrix. Following that, you could do something like this:

library(ggplot2)
library(network)

neto <- network(adj_mat, directed = FALSE, loops = FALSE, hyper = FALSE)
#assuming that it is undirected here and a simple network diagram.
#run ?network for more details in case your analysis is to be more complex

color_vec <- as.factor(color_vec)

col_pal <- c("blue", "red", "green") #assuming you only have 3 colors..

names(col_pal) <- levels(color_vec)
neto %v% "color" <- as.character(col_vec)

ggnet2(neto, node.size = .9, color = "color", palette = col_pal,
       edge.size = .05, edge.color = "white", mode = "adj")


how to filter out cell after doing QC I was wondering that the threshold that we set after doing QC in Seurat. there is an example

# GenePlot is typically used to visualize gene-gene relationships, but can be used for anything
# calculated by the object, i.e. columns in object@meta.data, PC scores etc.
# Since there is a rare subset of cells with an outlier level of high mitochondrial percentage
# and also low UMI content, we filter these as well
par(mfrow = c(1, 2))
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "percent.mito")
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "nGene")`

Then its filter cell command is under:

# We filter out cells that have unique gene counts over 2,500 or less than 200
# Note that low.thresholds and high.thresholds are used to define a 'gate'
# -Inf and Inf should be used if you don't want a lower or upper threshold.
pbmc <- FilterCells(object = pbmc, subset.names = c("nGene", "percent.mito"), low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))`

So I want to know whether the threshold setting is depending on the picture I show?

And the below is the Geneplot I have got base on my data. How can I set my thresholds?

Yes, you end up placing filters where things become inconsistent or similar to background. In the example, which appears to be from the Seurat tutorial, there's a downward plunge in nGenes vs. nUMI at around 200 genes and few believable cells with more than 2000-2500 genes, so the tutorials sets filters around there. The 0.05 mitochondrial contamination is a more arbitrary round number and is less supported by the plots (one would use more like 0.06 or 0.07 if one based this solely on the plots).


Entrez or Ensembl gene IDs? I'm using several datasets that are encoded using either Entrez or Ensembl IDs to specify genes, and need to decide on which to standardise on.

Are there any major reasons to use one over the other?

Which tools and conversion tables are best to use for this?

I think there is no general answer to your question.

Which type of ID you will use depends on how practical it will be for your downstream processing. This depends on what other source of information you have, how much conversion work you will need to do in order to match your data to these information if you use one or the other.

I think the choice may emerge by trial and error, as you realize that a given choice is actually not convenient in your particular use case. Try to set up easy to use conversion procedures that you can easily insert in some steps of your workflow, so that you can easily switch from one choice to the other.

Concerning the tools to do the conversion, see the following question and its answers: Converting gene names from one public database format to another


How can I find the position of every mutation where the Allele Frequency is greater than X, in regards to a reference such as Hg19? I have a Human genome sequenced in a BAM file (along with other files with the indels, snps, cnvs).

I want to find every mutation with regards to the reference Human genome. However the majority of "variants" are noise/read errors. I am therefore required to filter the variants via a minimum threshold based on allele freq (eg. 0,5).

Using a reference (eg Hg19), how can I find the position of every mutation where the allele frequency is greater than the threshold?

PS. I am just starting, so I don't know which are the best (command line?) tools for the job nor how to build a pipeline for this.

Help would be greatly appreciated

Hello and welcome to the world of bioinformatics,

BAM files are compressed formats of SAM and you can read more about it here.

To answer your question, I presume that your BAM is generated from a Whole Genome library preparation. Calling mutations on a Whole Genome takes a longer time than a Whole Exome since you are looking at every position instead of the exons of your genome.

What you are looking for is a special format of files called Variant Caller Format (VCF) generated by Variant Callers tools listed below in comparison studies. There isn't a tool better than another, they all have their advantages and disadvantages. I suggest that you read about them as they also describe their pipeline approaches.

Comparison of three variant callers for human whole genome sequencing
Comparing the performance of selected variant callers using synthetic data and genome segmentation
Performance assessment of variant calling pipelines using human whole exome sequencing and simulated data
Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers
Variant Callers for Next-Generation Sequencing Data: A Comparison Study

Once the VCF is generated, you will also have a column named AF which corresponds to the allele frequency also computed when using these tools. There are tools that allow you to filter VCF files to extract the needed information with a wished criteria like SnpSift,  vcftools, VariantFiltration. So, I suggest reading the above articles as a start.

As for building pipelines, it depends on what language you feel at ease at. However, if I could suggest, Nextflow is a good start. You can also use Snakemake if you like Python. Or maybe simply bash ? In any case, there are many pipelines already made to help you through your needs.

(If this helps, don't forget to either accept the answer, vote the answer or both).


Which one is a more convenient assembly? I have developed a software for de novo genome assembly. Its performance varies gradually according to how much data you employ. At initial stages it often produces contigs that look like that when aligned to a reference genome:

Initial-stage s-aligner assembly of SRR8357422
Genome fraction:98,8% -
Largest Alignment: 61.544 -
NGA50: 59.409 -
LGA50: 2
If I add more data, the contigs end up being larger, but that obviously takes more time. My question is, how does this result compare to, for example, one like that?

SPAdes assembly of SRR8357422
Genome fraction:93,46% -
Largest Alignment: 128.391 -
NGA50: 128.391 -
LGA50: 1
Note how you can construct a contig covering almost the full reference genome assembling green contigs in the first case. Therefore, a blast search of genes will likely let you find more genes given that it covers more genome fraction.
Apart from determining what option is more convenient, I would also like to know how much more convenient is to have both results. Or in which cases it's more convenient.
Last (but not least!), is there a way to generate automatically larger contigs from these overlapping contigs in assembly 1? I have tried using miniasm with no success. Even minimap2 seems not to find overlaps well enough. I am developing my own algorithm also for that but would like to know if there is another way.
[EDITED on 17/12 to add another example with real data for both assemblies]

I asked that on Reddit. I didn't get an answer to all my questions but I got an answer for how to assemble the overlapping contigs.
If I use Flye I can assemble these and get larger contigs, better NG50, and higher genome fraction.

s-aligner + Flye assembly of SRR8357422 Genome fraction:95,95% - Largest Alignment: 134.103 - NGA50: 134.103 - LGA50: 1
Given that the data in the first assembly apart from being more useful itself in my opinion, is objectively more suitable to generate a better result in all senses with a simple additional step that takes a few seconds, my only possible conclusion is that it is overall a better result and therefore this sort of result will often be preferable.
Of course, there are other aspects to consider, like error rate, that are not considered here as it wasn't part of the question (mainly because I think Quast is not adequate for calculating it). But BTW, mismatches per 100kbp are also lower.
[EDIT] I want to add that this is not granted that in all cases Flye will generate a clearly better assembly.


How to download the whole BLAST nt database into a specific folder? I have successfully downloaded the whole nt BLAST database into the current folder using:

wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"

However, I would like to download it in a specific folder, e.g. output/

I have tried:

wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz"  -O output/

but it does not work. Nothing happens. My 'wget-log' is empty.

I could move to that folder and download it there but I would like to remain in the same working directory.

C.

-O is to specify the name of the output file. You want the -P option.

So, from your example

wget -b "ftp://ftp.ncbi.nlm.nih.gov/blast/db/nt.??.tar.gz" -P output/


Differential Expression I have genes differentially expressed between two groups (case and control). I would like to annotate them by classifying them according to their biological functions.

I work on parasites of Plasmodium Falciparum type Dd2 which are divided into two groups (case and control). I carried out a Differential expression analysis with cuffdiff and CummeRbund then I annotated my genes on plasmodb.org.

Example of id genes: pfDd2_14005040

Now I would like to do a pathway analysis to see which biological pathways these genes are involved in.

StupidWolf,
I resumed my annotation with the 3D7 strain of plasmoduim Falciparum then I created my geneList as you indicated but when I make the dotplot nothing is displayed

setwd("C:/cygwin64/home/DIANGO/EXCELL")
library(org.Pf.plasmo.db)
library(clusterProfiler)

gene <- c ("PF3D7_0936800","PF3D7_1478900","PF3D7_1009700","PF3D7_0508500")
geneList = keys(org.Pf.plasmo.db)
enr <- enrichGO(gene,universe=geneList,OrgDb='org.Pf.plasmo.db',keyType ="SYMBOL",ont="BP")
head(enr@result,readtable=TRUE)
dotplot(enr)

The output gives a white rectangle

                   ID                           Description GeneRatio  BgRatio     pvalue
GO:0006323 GO:0006323                         DNA packaging       1/2  21/2615 0.01599974
GO:0071103 GO:0071103               DNA conformation change       1/2  30/2615 0.02281728
GO:0051276 GO:0051276               chromosome organization       1/2  58/2615 0.04387582
GO:0006996 GO:0006996                organelle organization       1/2 144/2615 0.10712138
GO:0016043 GO:0016043       cellular component organization       1/2 203/2615 0.14925925
GO:0006464 GO:0006464 cellular protein modification process       1/2 265/2615 0.19244223
            p.adjust     qvalue        geneID Count
GO:0006323 0.1026777 0.07205456 PF3D7_0508500     1
GO:0071103 0.1026777 0.07205456 PF3D7_0508500     1
GO:0051276 0.1316275 0.09237015 PF3D7_0508500     1
GO:0006996 0.2203602 0.15463872 PF3D7_0508500     1
GO:0016043 0.2203602 0.15463872 PF3D7_0508500     1
GO:0006464 0.2203602 0.15463872 PF3D7_1009700     1
> dotplot(enr)

Once again I need your help.


How to install/import Biopython into Python 3.8/ PyCharm IDE I'm new to Python/Biopython programming. I'd like to install Biopython into my PyCharm IDE but I've tried all instructions on this, but it's not confirming it. I have the available biopython file already extracted. There's a SyntaxError for 'install' in:
pip install biopython
What's the proper way of installing and importing biopython? And how does one confirm the installation.

Install Anaconda3 the instructions are here. Then making sure you are connected to the internet type
conda install -c anaconda biopython

Once Anaconda is installed don't attempt to use pip, if you do then you need to set an 'environment' in Anaconda and work from that environment.
I suspect you used pip or pip3 via
pip install biopython


Changing the legend names ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) + geom_boxplot() +
        scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3")) + theme_classic() +
        labs(title="Comparison",x="Norm", y = "counts") +
        scale_fill_continuous(labels = paste("condition1", "condition2", "condition3"))

I am trying to change the labels of the legend. Currently, it shows the labels that are present in my dataset. However, I want to change it to "condition1", "condition2", "condition3". Above is the code I am trying. It gives me following error.

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale. Error: Discrete value supplied to continuous scale In addition: Warning message: Removed 18 rows containing non-finite values (stat_boxplot).

Could someone help me here?

You don't have any continuous fill scales, since you're filling according to Condition, which is a factor. This is leading to the Error: Discrete value supplied to continuous scale error. You just need to combine your previous scale_fill_manual(values=...) command with the labels= part of your incorrect scale_fill_continuous() command:

ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) +
    geom_boxplot() +
    scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3"),
                      labels = paste("condition1", "condition2", "condition3")) +
    theme_classic() +
    labs(title="Comparison",x="Norm", y = "counts")

Something along those lines should work.


How to get the count of each kmer past 255 using khmer I have a Fastq file and I want to get the exact count of each possible kmer from this file.

On a previous post called How to use khmer to count k-mers? Daniel Standage proposed a custom script based on khmer methods that you can see down below :

>>> outfile = open('outfile.txt', 'w')
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz'):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

The code works if I switch counttable.get(kmer) with counts.get(kmer) but the main problem is that when I try on a customized fastq file with 300 occurrences of one kmer, the count will still be stuck at 255 and can't go further, and for some reason I can't pass through this limit.

Is this coming from the minimum table size ?

Does anyone know if there is a parameter that can work on a custom script to disable this limit ?

255 is the default maximum size of a Counttable in khmer. You want to do the following:

import khmer
counts = khmer.Counttable(31, 1e7, 1)  # Feel free to change the parameters
counts.set_use_bigcount(True)

The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (at some point you have to increase the size of the number used in the hash table).


Understanding PCHeatmap outputs I am currently trying to understand the purpose of these PCHeatmaps - part of the seurat package in R:

All the online documentation I have searched for has only highlighted how it is used to observe heterogeneity in the data.. This much I understand.

My questions are:

Are the genes in the rows of each PC plot just the most variable x
number of genes defined in FindVariableGenes?
If not, then how are each of these genes 'chosen' for use in the PCHeatmap
visualization?
What would be the different conclusions/deductions drawn from the PC1 heatmap and PC2 heatmap with regards to heterogeneity for example?
And finally, how can these conclusions help in the analysis of scRNASeq data?

The PCHeatmap function (renamed DimHeatmap in Seurat v3) can be used to help determine the number of principal components to use in downstream analysis, as well as to visualize the top genes contributing to each PC.

Both cells and genes are ordered by their PC scores, and by default the 15 genes with the highest and 15 genes with the lowest PC loadings are displayed for each PC.

In your case with PC1 and PC2, my conclusion would be that PC1 and PC2 both strongly separate different populations of cells. For PC1 there is a smaller population with very strong separation, and PC2 splits the cells roughly in half.

You can also set the cells.use parameter (renamed cells in Seurat v3) to only show the top n cells for each PC which can help with visualization.


Produce a single sequential FASTA sequence out of BAM I'm having problems properly looking for a solution because I'm a layman in Bioinformatics not familiar with the terminology. I'm hoping you can nudge me in the right direction, please! Thank you very much.

What I want in the end

A FASTA file containing a single sequential sequence (chromosome 1-22..) of the human genome.

What I have

23gb of "clean_data" which are .fq.gz files.
21gb of "result_alignment" which is .bam+.bai
"result_variation" indel.vcf.gz and snp.vcf.gz

What I tried

I used samtools fasta reads.bam > reads.fasta which results in many lines like so:

>V100005418L4C001R001322578/1
AAATAGCCCACAC...
>V100005418L4C001R006084507/1
GCCATAAAGCCTA...
...

But these include many many duplicate lines and the reference genome from the human genome project GRCh38 is only 3GB in size, so I thought that this FASTA sequence needs to be properly "aligned" somehow?

Am I thinking correctly here?

I'm really confused with all of the terminology and don't know exactly what I need to look for. Could someone kindly put the process into simpler words for me?

If you just want the human chromosomes in a FASTA format, why not downloading directly the individual chromosome files (chr*.fa.gz)?

If you don't rely on the "official" releases (for whatever reason), then what you need is first to assemble your reads into genomes, possibly (but not surely) going down to a chromosome-level resolution. Different approaches for assembly are possible/available, mainly depending on (1) what do you biologically need to answer and (2) the type of sequencing data you have - we'll need more details to be able to help you in this direction.

The samtools fasta utility will "just" convert your reads (likely given from a sequencing facility) from an alignment (BAM) format into a reads (FASTA) format, there is no concept of chromosomes here. You won't get assembled chromosomes this way, just extract the sequence in FASTA format from your alignment file.


How to remove frame on RDKit figures? RDKit can plot molecules, thought the structures are surrounded by an ugly frame with ticks. How can I turn that off and plot the molecules without that frame?

Note, the %pylab inline portion of this code is from the Jupyter environment.

%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')

It looks like rdkit and rdkit's Chem just use the matplotlib. So, if you import matplotlib you can edit your figure like any other matplotlib figure.

I produced this from the command line, so it may work differently in the Jupyter environment.

It looks like Chem doesn't properly scale the drawing based on how large the chemical is, so I added plt.xlim and plt.ylim to make the molecule larger in the frame.

from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
plt.title('Test')
plt.axis('off')
coord_min = 0.35
coord_max = 0.65
plt.xlim((coord_min, coord_max))
plt.ylim((coord_min, coord_max))
fig.savefig('mol.jpeg', bbox_inches='tight')


How can I get bedtools to tell me which genes are being expressed? I'm trying to align the Acinetobacter baumanii genome to a genome. I've already done the alignment, and I want to use bedtools to see which genes are being expressed exactly. When I try running the following command,
bedtools coverage -a outCSEQSAligned.out.bam -b GCF_000746645.1_ASM74664v1_genomic.bed

unfortunately, it doesn't tell me which genes are being expressed at each place but instead tells me that the sequence aligns to the single Acinetobacter chromosome which isn't too helpful:
NZ_CP009257.1   282134  282254  GWNJ-1012:340:GW200702637th:1:1105:15519:10019/1        255     +       282134  282254  0,0,0   1       120,    0,4       120     120     1.0000000

I'm assuming that 282134-282254 are the positions on the chromosome that correspond to that alignment. I suppose I could extract these positions and figure out which genes they correspond to, but I'm wondering if there's an easier way to use bedtools to figure out which genes are being expressed. For reference, outCSEQSAligned.out.bam is my bamfile which I aligned the Acinetobacter genome to, and GCF_000746645.1_ASM74664v1_genomic.bed is a bedfile corresponding to the Acinetobacter genome.
Thanks for the help!

Bedtools is not the best tool for this. You can use featureCounts, which gives you the number of read that overlap with features in a genome. In this case, your features are genes.
bedtools coverage is used to calculate the overlap between each read in BAM and each line in BED. It is not very straightforward to quantify gene expression. Plus bedtools does not account for exon in each gene and summarize them to gene level, like featureCounts does.


What is the recognized way to reference accession data (.sra files) from NCBI as a URI? I would like to assign a URL/URI to accession runs (i.e. SRR IDs) from the NCBI short read archive.

In program listings, I used to list them as ftp:// URLs, because they were available on their servers as such, e.g. SRR3418242:

ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR341/SRR3418242/SRR3418242.sra

But this approach is deprecated, plus, not all runs are available this way anymore. The now recommended approach is to use sra-tools to download the data. This tool takes the SRR Ids directly, i.e. prefetch "SRR3418242", and resolves them to one of many possible locations (home directory, http, etc).

Is there an official URI scheme that can be used to reference accession runs from the short-read-archive ? I would like to disambiguate them from other types of files and URLs (S3 keys, http, file). For file manifests, spreadsheets, and publications, is there something similar to srr://SRR3418242?

The nomenclature is no longer being followed as NCBI is moving all the files to S3/Google cloud. pysradb allows you to fetch the metadata and get the URLs for each individual SRR. See examples in this notebook.

You can also download an entire project by doing:

pysradb download -p <SRP>


What are the coordinates of the extended major histocompatibility complex (xMHC) in GRCh38? Is there an accepted boundary of the extended major histocompatibility complex (xMHC) in human? Specifically, I am interested in the boundary coordinates for Genome Reference Consortium Human Build 38 (GRCh38 / hg38).
I'd like to use these coordinates to query MyGene.info to retrieve a list of genes in the xMHC. For example, King et al 2019 used the following definition:

we excluded drugs with nonhuman or xMHC gene targets. We considered xMHC to include HIST1H2AA and KIFC1 and all genes between them (Chromosome 6 25.7 Mb-33.4 Mb).

But I was curious whether there was an authoritative definition of the xMHC region in humans that would specify down the the base pair rather than just megabase pair.

Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13
which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.
The "extended MHC" defined here [1] is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The xMHC is defined in [1] as the regions between HIST1H2AA and RPL12P1 genes, which would give you chr6:25726063-33400644 in Grch38. This is very similar to your HIST1H2AA - KIFC1 definition (chr6:25726063-33410226). Of course if you're trying to reproduce a result, probably better to use the original paper's definition.
[1] https://www.nature.com/articles/nrg1489


Uploading files from local computer or web in rstudio cloud I am using RStudio Cloud for running velocyto R package. So I have tried to import velocyto file from http://pklab.med.harvard.edu/velocyto/mouseBM/SCG71.loom

I tried to upload that like so

But saying

> SCG71 <- read.table("/cloud/project/SCG71.loom", header=TRUE, quote="\"")
Error in make.names(col.names, unique = TRUE) :
  invalid multibyte string at '<89>HDF'
In addition: There were 19 warnings (use warnings() to see them)

Can you help me in very initial of velocyto R?

This seems to be an encoding problem (certain characters not being recognized), see this, this and also this. Probably the file "/cloud/project/SCG71.loom" is encoded with a format that is not included in your "locale".

Check the output of Sys.getlocale() and see if your input file's encoding is listed. If not, you will need to try something like Sys.setlocale("LC_CTYPE", "en_US.UTF-8") or Sys.setlocale("LC_ALL", "English") (see the third link for details).


can I download DESeq2 in R 3.6.3 in Linux MInt? I am trying to download DESeq2 in R 3.6.3. Is it possible?
This is printed when I am trying: Warning in install.packages : package ‘DESeq2’ is not available (for R version 3.6.3)
Thank you!

Bioconductor is a package management tool for R that allows you to install the appropriate package versions for your R version.
For R 3.6.3 you'll need Bioconductor 3.10:
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(version = "3.10")

Then you install  DESeq2:
BiocManager::install("DESeq2")


Best practice to perform a functional analysis on a set of VCFs from a cohort of patients I am confused about how to perform a functional analysis after a variant calling on a cohort of patients. I have anotated all vcfs and filter them in order to get only those deleterious variants.

I was thinking on merging those vcf:
 I tried to merge them with bcftools merge and also using VariantAnnotation and transforming each sample vcf into a dataframe containing the variant lines to then row bind all those dataframes together.

As an example, Ill show some fields of teh first row of a dataframe:

CHROM   START REF ALT FILTER DP  GENEID  SIFT
  17     1000   T   A   PASS   62   TP53    0

Then, I would extract the genes affected by those variants and perform a gene set enrichment analysis, a pathway analysis similar to those performed after an RNA-seq analysis (with enrichR, for example). But I am confused about the viability of this steps. For example, I would like that, during the pathways analysis, the program I use not just reads the list of mutated genes in my cohort, but also takes into account that, if a gene is mutated in more than one sample, it should wheight more in the analysis than a gene which appears in just one sample.

As a simple approach you may want to try the "fuzzy" list implemented in EnrichR itself. There is some information about this in the Help Center under the "What is a gene set" heading.

This way you can code the number of data sets a gene appears as a value between 0 and 1. In the Help Center this is used to convert fold changes or p-values to a membership score but I don't see why this would not also work for the number of data sets a SNP is observed in.

Since it is relatively easy to implement this can give you an first impression of what to expect. It then may be combined with analyzing each data set independently and then comparing the overlap.


Help with setting DimPlot UMAP output into a 2x3 grid in Seurat On my merged seurat object of 6 samples, when I use the split.by function in tandem with the Dimplot/UMAP plot, all six samples are displayed in series along a commonly labeled 'UMAP_1' x-axis in an order of (these are arbitrary for simplicity, here) Z, Y, X, C, B, A. Instead, I would like to specify these plots to output as a 2 row x 3 column arrangement ordered (1st Row) A, B, C and (2nd Row) X, Y, Z.
I'm assuming I could use something like: group_levels <- c("A", "B", "C", "X", "Y", "Z") to define the order I'd like, but I'm lost on how to get this to arrange the UMAP plots into a 2x3 grid I'd like.
Appreciate any help!

how to get this to arrange the UMAP plots into a 2x3 grid

You can use the ncol parameter, which will set the number of columns in the "combined plot".
And for future reference, if you get the help page of the function you are using with ?DimPlot, you will see the .... This means that you can feed additional arguments to other functions that DimPlot() uses. Apparently,  CombinePlots() is used under the hood and has an ncol parameter.


Preparing PDB file for modeling with Swiss-Model I recently starting studying bioinformatics and have a question. I need to build a model of protein using Swiss-Model. First, I need to download a protein record from PDB (for example, 2HZI), then I must open PDB file using simple text editor. I know this usually requires biopython, but i don't use it yet. Now I'm trying to study how to build a model.

Swiss-model needs only one chain to build model, and I remove from file all extra information like AUTHOR,REMARK,HETATM etc. But Swiss-model still gives an error: "Sequence String contains illegal characters. Only letters and gaps allowed."

As explained in the website the data introduced must be:

  Format must be FASTA, Clustal, plain string, or a valid UniProtKB AC

So you need to paste just the protein sequence. Without seeing the data you paste I would guess that you don't use a valid plain string.

As Jonathan Moore explains below a FASTA format is:

>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL

(breaking the sequence in several lines is optional but usually done fore easy reading)


Retrieve RNA sequencing data for human p53 colon cancer cell lines I want to download the next generation RNA sequencing data in different cell lines. I want the data (alongwith normal) for p53 wild type, p53 knock-out and p53 null cell lines (human colon cancer). From where do I download this data? In fastq format if possible.

Currently, I am going through SRA, to look for the datasets, using their search engine. But the process is extremely tedious.

Also can someone please hint about how to download the data from SRA, in case that's the only way to get data. I'm getting very confused.

I suggest you to try difference resources, including ArrayExpress.
A simple search at ArrayExpress returns multiple experiments which might be of your interest. Just a pointer, ArrayExpress is also synced with GEO.


Downloading dataset from SRA (SOLiD Platform) Trying to download colorspace data from SRA, but getting an error

abi-dump -A SRR1657115.sra

abi-dump.2.9.6 err: item not found while constructing within virtual database module - the path 'SRR1657115.sra' cannot be opened as database or table.

You have to use the accession name without the .sra suffix.

abi-dump -A SRR1657115


Evidence of emergence of "genuine" novel protein interactions by aa mutation Anyone knows of an example of a protein that, without coming from a recent duplication event, underwent a mutation(s) that caused it to have a novel interaction with a new ligand, substrate, other protein or molecule?
By novel I imply that it doesn't "go back" to an ancestral interaction or doesn’t develops an interaction that it's already performed by a close relative.
I assume that these mutations are extremely rare, and that a duplicate is almost always necessary to have them fixed in a lineage.

HA Sialic acid receptor binding in influenza virus determining host specificity, viz. Bird flu vs human flu.


semantic similarity measurement for cell line ontologies I have a set of cell line pairs and I want to know to what extent the pairs are similar based on their ontologies.
The problem I have is that I have found a Python library called Fastsemsim, but it has many similarity measurements. Or to be exact, there are actually many ontology similarity measurements. And I want to know which similarity measurement is best suited for my case (cell line ontologies similarities). There might be a couple of good ones that are good for cells ontologies.
These are a set of similarity measurements that are available in Fastsemsim.

Semantic Similarity

SemSim.TermSemSim
SemSim.ObjSemSim
SemSim.ObjSetSemSim
SemSim.SetSemSim
Specific Semantic Similarity measures
SemSim.MixSemSim
SemSim.avgSemSim
SemSim.maxSemSim
SemSim.BMASemSim
SemSim.SemSimUtils

Source:
https://fastsemsim.readthedocs.io/en/latest/library.html

About the metric you use, there are several reviews on the topic, if I recall well I think that the information content based are the best ones, of those, I only recognized the JiangConrathSemSim listed on the webpage.

However, if you want to use a semantic similarity measure I would use the BMA method, which is more comprehensive than other methods. But the BMA method is only a way to aggregate several metrics of ontologies into a single numeric value. I am not confident that FastSemSim implements correctly those methods.

Note, that according to the license, you should contact the authors before using it:

  You must contact us before using it or its results or any application/tool based on it in any published work.

I am not sure that the cell lines fit in a ontolgogy. Because ontologies are defined as directed acyclic graph (DAG) and there are some cycles that happens in cells like the induced pluripotent stem cells show. However they might be right if you just compare the cell lines and not the state of the cells on those lines (use it only if you want to compare the lineage of your cells lines).


NCBI SRA database sample: control vs test I was trying to download some data from NCBI SRA (SRA059451). There are 27 samples available for SRA059451. But i am unable to understand which samples are 'control' and 'test' samples. Please help me regarding how to differentiate the samples.

You can download the information using this link:

I did it in R:

tab = read.delim("https://www.ebi.ac.uk/ena/data/warehouse/filereport?accession=SRA059451&result=read_run&fields=run_accession,experiment_title,sra_ftp,sample_alias")

head(tab)
  run_accession
1     SRR585697
2     SRR585700
3     SRR585726
4     SRR585727
5     SRR585728
6     SRR585729
                                                                experiment_title
1 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 171
2 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 238
3 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 241
4 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 242
5 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 243
6 Illumina HiSeq 2000 sequencing; Symptomatic atherosclerosis gut metagenome 245
                                      sra_ftp
1 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585697
2 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585700
3 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585726
4 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585727
5 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585728
6 ftp.sra.ebi.ac.uk/vol1/srr/SRR585/SRR585729
                                    sample_alias
1     Symptomatic atherosclerosis gut metagenome
2 Symptomatic atherosclerosis gut metagenome 238
3 Symptomatic atherosclerosis gut metagenome 241
4 Symptomatic atherosclerosis gut metagenome 242
5 Symptomatic atherosclerosis gut metagenome 243
6 Symptomatic atherosclerosis gut metagenome 245

Normally the information will be here, but I cannot find what does metagenome 238 etc refers to. No table in the paper that can link sample to phenotype..

Most likely you need to contact the authors. This list might be useful though.


Drug with as many targets as possible through a system biology graph? I have 4 pandas data frames containing the following data:

affinity_score1: affinity between each couple (Drug-Target) measured using  tool1
affinity_score2: affinity between each couple (Drug-Target) measured using tool2

The two tools (tool1, tool2) give measure in different scales and are not 100% accurate.

Drug_similarity: a similarity matrix between drugs
Target_similarity: a similarity matrix between targets.

Here also the tool used to measure the similarity is not 100% accurate.

import pandas as pd

affinity_score1 = [{'Drug0': -24.5, 'Drug1': -17.1, 'Drug2': -13.1},
           {'Drug0': -20.0, 'Drug1': -21.6, 'Drug2': -18.1},
           {'Drug0': -19.0, 'Drug1': -15.4, 'Drug2': -20.0}]

affinity_score2 = [{'Drug0': 1.5, 'Drug1': 0.8, 'Drug2': 0.4},
           {'Drug0': 1.2, 'Drug1': 1.3, 'Drug2': 0.9},
           {'Drug0': 0.8, 'Drug1': 0.7, 'Drug2': 0.9}]
affinity_score1_df = pd.DataFrame(affinity_score1).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
affinity_score2_df = pd.DataFrame(affinity_score2).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})
Drug_similarity = [{'Drug0': 1, 'Drug1': 0.8, 'Drug2': 0.2},
           {'Drug0': 0.8, 'Drug1': 1, 'Drug2': 0.9},
           {'Drug0': 0.2, 'Drug1': 0.9, 'Drug2': 1}]
Target_similarity = [{'TargetX': 1, 'TargetY': 0.8, 'TargetZ': 0.2},
           {'TargetX': 0.8, 'TargetY': 1, 'TargetZ': 0.9},
           {'TargetX': 0.2, 'TargetY': 0.9, 'TargetZ': 1}]
Drug_similarity = pd.DataFrame(Drug_similarity).rename(index={0:'Drug0',1:'Drug1', 2:'Drug2'})
Target_similarity = pd.DataFrame(Target_similarity).rename(index={0:'TargetX',1:'TargetY', 2:'TargetZ'})

I am trying to find the Drug who can have as many Targets as possible.

So far, a simple approach could be to select the Drug having the highest average score in both affinity scoring tools.

This approach does not use the information in the similarity matrices.
Given that if a Drug1 is successful with TargetX, and TargetX is similar to  TargetY then Drug1 can also easily be successful with TargetY. Considering the inaccuracy of the scoring tools, I think including this information could improve the accuracy of the answer, eliminating some false positives.
Also, if a Drug from a cluster of similar drugs (information that can be derived from the drugs similarity matrix) is the only one predicted to have good affinities with some targets, then this is more likely to be a false positive.
But I don't know how to efficiently integrate all the information. Scipy may help with clustering the two similarity matrices.

Another way of analysis is to set affinity and similarity thresholds, in which case each data frame could be transformed into a network:

Each node is a drug or a target.
Two types of edges: similarity (between drugs or targets) and affinity (drug-target)
The connected nodes have passed the thresholds.
The edges have weights: affinity or similarity score.

The problem then becomes finding the most connected node (a drug node) with the highest weight values (scores on tool1 and tool2). Connected nodes through similarity may then be used to validate the number of targets for each drug.
It seems a more complicated approach but more rigorous with most likely higher positive predictive value.

Any python package to implement such an approach?
And also a visualization tool.

The original data has 38 Targets with 806 drugs.

I don't know any python (or R package) to do this), but this is how I would do it based on these matrices:

Do a multidimensional scaling (MDS) plot of the target_similarity matrix.  Explore the visualization:
Find which drug targets the drug closest to the 0,0 point in the 2D space

However this doesn't say anything about the affinity between the drug and the target.
If you wanted to add this then you could weight how much close is the distance between the drug and the target vs how close are they to the center of the MDS (a simple division could work if the distance is scaled between 0 and 1 and averaged).

Note that depending on the original information used to calculate the similarities and the affinities it could be done in other ways


Interpreting 0x200 flag in bwa-mem alignments I am looking at the bwamem.h code in http://github.com/lh3/bwa and found that BWA-MEM will give flag 0x200 to what it calls MEM_F_SOFTCLIP:

#define MEM_F_PE        0x2
#define MEM_F_NOPAIRING 0x4
#define MEM_F_ALL       0x8
#define MEM_F_NO_MULTI  0x10
#define MEM_F_NO_RESCUE 0x20
#define MEM_F_REF_HDR   0x100
#define MEM_F_SOFTCLIP  0x200
#define MEM_F_SMARTPE   0x400
#define MEM_F_PRIMARY5  0x800
#define MEM_F_KEEP_SUPP_MAPQ 0x1000

But looking at the explanation in this website: https://broadinstitute.github.io/picard/explain-flags.html

It labels 0x200 as read fails platform/vendor quality checks.

If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC? Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?

The SAM/BAM specification labels 0x200 the same way as the Picard flag explaining website:

The question remains as to why BWA-MEM uses the 0x200 flag for softclipping notation rather than Fail QC.

bwa mem is using exactly the same meaning of the 0x200 flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for primary and secondary alignments (hard-clipping the latter).

For what it's worth, bwa mem (at least the version I just tried) doesn't set the QC failed flag. This behavour is actually pretty common among aligners, since Illumina has changed how they represent QC failure at least once.


Finding data to validate results My group has a complex datset and I would like to validate the results and check the methods I use in a similar dataset (better if it has been already studied).
Characteristics of my dataset:

It is related to the inflammatory bowel disease (IBD)
More specifically to Crohn's disease.

It has RNA-seq
I have explored the GEO and there are few projects (10) with RNA-seq on the tissue of my disease, and I don't know how to find the 16S-seq of those projects if it exists (see point below).

It has 16S-seq (and not only from the stools)
The ENA seems to have 16S sequencing but most of them are from stools not from biopsies and still, I don't know how to find the RNA-seq of those projects if it exists (see point above)

It has both RNA-seq and 16S-seq from the same individuals at the same time point.
The NIH Human Microbiome Project for my disease has only one individual with that both 16S-seq and RNa-seq of the same patient at the same time.

How can I find a dataset similar to mine?

The snarky answer is that it's better to find a matching dataset at the experimental design phase (i.e. before these studies are carried out), so that the bioinformatics component on the other side of things ends up being a lot easier and more robust. Although it can be quite time consuming to do a database search at the start, it's even harder to find matching projects if the legwork hasn't been done during the design phase. Researchers should be thinking about how hypotheses will be validated before the experiments are carried out.

Every additional characteristic that is chosen for comparison adds complications with regards to matching data, and reduces the applicability of results when there is a mismatch. You would have to be somewhat lucky to be able to find a IBD-Crohn's dataset with paired RNA-seq data and 16s-seq (presumably short-read; is 454 good enough?) from a specific non-stool tissue.

Apart from the places you've already looked, there's dbGaP and SRA. Given the number of characteristics that have been mentioned, I expect you'll have to make some compromises when comparing other datasets.


Visualizing cell growth I model the following events: a rod-shaped cell with center $x$ grows symmetrically until it reaches a maximum length/age, it divides into two identical cells, and the process continues. At division, the direction angle $\varphi$ of the daughter cells changes randomly.

So, the growth rate is $l=\frac{l_{max}}{2}+\frac{a}{a_{max}}\frac{l_{max}}{2}$

Is there any way of naively visualizing this? I found different cell modellers on GitHub but they were very complicated biologically speaking. I have some basic Python knowledge, but I have never done simulations.

I would recommend looking into Python's plotting utilities, for example here.

I would recommend specifically looking at the examples that plot mathematical functions.

The basic idea is that you create a function that computes your equation based on the different parameters as inputs ($a, a_{max}, l_{max}$) and then you can plug in a range of values for each parameter and plot e.g. with $a_{max}$ on the x axis, then different lines showing a range of values for $l_{max}$, or whatever makes sense for you.

e.g.

import numpy as np
import matplotlib.pyplot as plt

def fn(a, amax, lmax):
  return lmax/2 + (a/amax) * (lmax/2)

plt.figure()
amax = np.arange(0.01, 1, .01)
plt.plot(amax, fn(.005, amax, .1), "ro") ## red points
plt.plot(amax, fn(.005, amax, .2), "go") ## green point
plt.show()

I have no idea if these parameter ranges are meaningful for you, you will have to figure that out for yourself.


Is there a tool for synonymous recoding of DNA sequences? I've got a DNA sequence that I'd like to make synonymous mutations throughout, thereby preserving the amino acid sequence. Does anyone know of a tool to achieve this

There seems to be a bioconductor package named synmut which does just that.
Apparently, it can take codon usage into account when generating synonymous mutations.

Copy pasting an example that matches your use case from the official doc:

library("SynMut")
filepath.fasta <- system.file("extdata", "example.fasta", package = "SynMut")
filepath.csv <- system.file("extdata", "target_regions.csv", package = "SynMut")
region <- read.csv(filepath.csv)
rgd.seq <- input_seq(filepath.fasta, region)
mut.seq <- codon_random(rgd.seq, keep = TRUE)
Biostrings::writeXStringSet(get_dna(mut.seq), "rgd.fasta")


Running SnakeMake on cluster I do not understand how to specify correctly the parameters on a SLURM cluster for snakemake to use them. I tried submitting the following SLURM file, but it does not work this way and the number of cores used is only 1, not 20:

#!/bin/bash

#SBATCH -p standard
#SBATCH -A overall
#SBATCH --time=12:00:00
#SBATCH --output=snakemake%A.out
#SBATCH --error=snakemake%A.err
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=20
#SBATCH --mem=120000

snakemake

Then, I tried following snakemake tutorial. And created cluster.json based on the SLURM parameters that I need:

{
    "__default__" :
    {
        "A" : "overall",
        "time" : "24:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 20,
        "p" : "standard",
        "mem": 120000,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

And ran snakemake inside a newly created snakemake.sh script:

#!/bin/bash

snakemake -j 999 --cluster-config cluster.json --cluster "sbatch -A
{cluster.A} -p {cluster.p} \
-t {cluster.time} --output {cluster.output} --error {cluster.error} --
nodes {cluster.nodes} \
--ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem
{cluster.mem}"

And it is giving me an error now:

  sbatch: error: Unable to open file
  /bin/sh: line 1: -t: command not found
  Error submitting jobscript (exit code 127):

I am now completely lost at what I shoould actually do. I would prefer plain regular .slurm file submission, but how to make snakemake use them? Any suggestions would be greatly appreciated.

I removed \ - line separators in the snakemake.sh script:

#!/bin/bash

snakemake -j 10 --cluster-config cluster.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --nodes {cluster.nodes} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

And it started to run. It is not convenient though for me. I would rather prefer submitting just one job using .slurm file passing all of the parameters from #SBATCH. Is it possible?

You can omit --nodes, you need the following:

#!/bin/bash
#SBATCH --ntasks-per-node=1
#SBATCH -c threads_from_snakemake
#SBATCH -p some_partition

Commands go here

For slurm, you may want to modify my SlurmEasy script. We use that with snakemake all the time, usually of the form:

snakemake -c "SlurmEasy -t {threads} -n {rule} --mem-per-cpu {cluster.memory}" ...

You'll also want units with your memory requests, like 12G. Also, in general it's best not to submit snakemake itself to the cluster, rather run it on an interactive node and have it submit jobs on your behalf (if the admins complain point out that snakemake is barely using any resources on the head node).


Pipeline (Pypiper) returns error, but individual command works fine I am running a script file on a cluster with slurm through sbatch. Locally it runs perfectly well on 2 different linux versions (linux mint, ubuntu 12.04), however, it returns -9 error code while running on a cluster which is CentOS 7.3. When we are using Pypiper and run the whole .py script file, it fails when it starts running trimmomatic program

http://www.usadellab.org/cms/?page=trimmomatic

However, when running just trimmomatic program step by itself it does not fail and finishes successfully. We have a suspicion that something may be wrong with how Pypiper runs on the cluster, however, feeling rather stuck.
Any suggestions would be greatly appreciated.

Here is the error that is printed out:

ILLUMINACLIP: Using 2 prefix pairs, 22 forward/reverse sequences, 0
forward only sequences, 0 reverse only sequences
</pre>
Process 109391 returned: (-9). Elapsed: 0:00:51. Peak memory: (Process:
5.211GB; Pipeline: 0.0GB)

### Pipeline failed at:  (12-14 15:50:47) elapsed:0:23:09 _TIME_

Total time: 0:23:09

Change status from running to failed
Traceback (most recent call last):
  File "...../salmon_rna_seq.py", line 344, in <module>
sys.exit(main())
File "...../salmon_rna_seq.py", line 304, in main
pipeline.process()
File "...../salmon_rna_seq.py", line 277, in process
self.trim_reads()
File "...../salmon_rna_seq.py", line 195, in trim_reads
self.pipe_manager.run(cmd, self.sample_files['r1']['fastq']['trimmed_paired'])
File "...../pypiper.py", line 508, in run
process_return_code, local_maxmem = self.callprint(cmd, shell, nofail, container)  # Run command
File "...../pypiper.py", line 691, in callprint
self.fail_pipeline(e)
File "...../pypiper.py", line 1111, in fail_pipeline
raise e
OSError: Subprocess returned nonzero result.

Pypiper terminating spawned child process 108475...
child process terminated
slurmstepd: error: Exceeded step memory limit at some point.

The last line may be informative. Although, we are still not sure what we could possibly do.

Slurm is killing the process because it's using more memory per core than the default allocation (that's what slurmstepd: error: Exceeded step memory limit at some point means). I've never used pyPiper (like seemingly most people, I use snakemake), but presuming that you have it running either srun or sbatch at some point then specify something like --mem-per-cpu 10G.


bedtools intersect on very large .bam file - -sorted confusion I need to do a bedtools intersect operation on a very large .bam file.  When I use the standard bedtools intersect operation, the process consumes all the memory on my system.

From the bedtools documentation page, I find the -sorted option to bedtools intersect, which sounds like it is the perfect solution to the problem, but the documentation seems incomplete.  From the documentation, I need to ensure that the .bed file is sorted and I can use either bedtools sort or the Linux sort command to accomplish this:

bedtools sort -chrThenSizeA -i input.bed > sorted.bed

sort -k1,1 -k2,2n input.bed > sorted.bed

If I then execute the intersect operation on the .bam file, I get an error.

bedtools intersect -sorted -c -a sorted.bed -b input.bam > intersect_linuxsorted.out
ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files.

Tracking this down I find that samtools sort, which was executed on the .bam file, "sorts by the leftmost coordinate, or the read name if -n is used."  It would appear there is no way to put the .bam file in order by reference name, which is what bedtools instructs.

For a workaround, bedtools sort does have the -g option, which allows me to specify the ordering of the .bed file.  To get this ordering, I execute:

samtools view input.bam | cut -f3 | unique > order.txt
bedtools sort -g order.txt > alt_sorted.bed

This puts the .bed file into the same order as the .bam file, and the bedtools intersect -sorted command now works, but is it working correctly?  It seems that bedtools intersect -sorted wants the files in consistent ordering, which can be done, but is that consistent with the intersect algorithm itself?  From the documentation's instructions to sort by the first 2 fields, chromosome name and left-most coordinate, it seems inconsistent to order it as the .bam file is ordered, by read name and not reference name.

I'm not sure how to check this assumption.

EDIT:  Some extra details regarding sorting of the .bam file.

If the .bam file is sorted with samtools sort input.bam > sorted.bam I get the error, 'ERROR: chromomsome sort ordering for file input.bam is inconsistent with other files`

If the .bam file is sorted with samtools sort -n input.bam > sorted.bam I get this error, Error: Sorted input specified, but the file sorted.bam has the following out of order

It would seem that if I sort without -n, it thinks the chromosomes are out of order, and if I sort with -n option it thinks the coordinate order is wrong. The only fix I have is forcing the order of the .bed file with the -g option, but again, I have no idea if that is consistent with the underlying algorithm's assumptions for the -sorted option.

It sounds like bedtools is behaving properly. The bam file is sorted by read name only if you used the -n option in samtools sort. Otherwise, if you used samtools sort without the name-sort option, it is 'ordered' by reference name (which can be arbitrarily ordered, as I'll explain below) and then by position (which is numerically ordered).

I think you are getting confused about 'ordering' the reference names. There are multiple ways to sort names; Some people might want numeric sorting (chr1, chr2, chr3...) and others might want lexicographic sorting (chr1, chr10, chr11...). Because of this, many tools allow you to specify the way you want to order the reference names.

There isn't necessary only one correct way to sort the reference names; what matters is that the order is the same between the files you're trying to intersect. This is exactly the purpose of bedtools intersect -g.

If you want to check the order of reference names in your bam, you can easily do so by checking the header via samtools view -h input.bam. It will say something like this:

@HD     VN:1.3  SO:coordinate
@SQ     SN:chr1     LN:248956422
@SQ     SN:chr10    LN:133797422
@SQ     SN:chr11    LN:135086622
@SQ     SN:chr12    LN:133275309

The @SQ lines tell you how your bam is sorted. My suspicion is that this order is different from the Linux sort or bedtools sort default, which is why your bedtools intersect was failing.


Encoding a 15bp DNA sequence into a shorter than 15 digits number I have a series of 15bp DNA sequences from single-cell barcodes and I would like to encode them into a shorter than 15 digits number so that each unique 15bp DNA sequence will create a unique numeric value. What's the best way of achieving this?

What you are asking is called a "hash" (wiki).
In python, get your oligo as a string and use the standard library method hash(oligo), this is a number that gets used for a lot of operations. It is however always a 18+1 digit number in decimal, so a different solution is needed.
There are four nucleobases, so binary conversion from scratch is easy and best option and each solution will be unique and small. Each base takes up 2 bits and you concatenate these and convert them to hex or decimal —the length should be 10 decimal digits. One issue is when you have a oligo of a different length which start with adenosine or not.
One issue with this is the order. I am not sure where it comes from but the bases are nearly often expressed in the order ATGC. So...
mapping = {'A': '00', 'T': '01', 'G': '10', 'C': '11'}

It is satisfying how the first digit is the W/S, while the second is R/Y.
Now simply covert and add them compressed = functools.reduce(lambda acc,base : acc+mapping[base], oligo, '') and covert the string to decimal like so int(compressed, 2).


Melt p-values for CpG sites mapping to the same gene I have some data I am working with, and I am curious if I am able to combine p-values from a paired t-test for CpG sites in the genome using Fisher's Method to get one p-value for each unique gene. Linked here is the Wikipedia page for Fisher's Method. I understand that an assumption of the method used is that each individual p-value being combined must be independent. I'm relatively new to biostatistics, so I'm curious if using CpGs from the same gene would violate this assumption.

Methylation levels have high local correlation, so Fisher's method would be problematic. Having said that, you have no reason to use Fisher's method after a paired t-test. A paired t-test will give you a single p-value per gene, which is what you want. Do be sure to only include CpG with some minimal coverage in both group.


What unit is the length of a protein measured in? Sorry if this question is obvious - I'm not familiar with proteomics so I'm not sure of the standard units.
In this paper, it talks about using the length of a protein to calculate the NSAF value.  However, it doesn't specify the units.  Would it be peptides, amino acids or nucleic acids? Or something else?

As mentioned by others, it's the number of amino acids. Alternative measurements could be in Angstroms or similar, but that would require at least an approximate structure.


For what bioinformatics tasks is Biopython more adapted than Bioperl? Are there any advantages to learning Biopython instead of learning Bioperl?

Ideally, we would learn both, but someone starting out in bioinformatics may have to choose what to learn first depending on the kind of problems actually encountered.

Are there problems for which Biopython is better than Bioperl (or vice-versa)?

Regading the perl vs python discussion, there is no final answer which language is better, but I have some advice for you:

Learn the language your colleagues or your advisor use. This way you are able to discuss your code with them and also get help if you run into problems.


How can I run the fdnadist command from the EMBOSS package? My operating system is Mac OS X.
I would like to compute DNA distances between different sequences. For this, the dnadist package from the Phylip software would be good but hard to automate. EMBOSS seems to have an fdnadist command for this. I have installed emboss version 6.6.0 from the bioconda package but if I type fdnadist, the command is not found. How can I use dnadist or fdnadist on the command line?

To find emboss packages type wossname, then hit enter twice.  You will see all the programs.  There is one called distmat.  Then do tfm distmat, to get the fkn manual for distmat.


nextflow `fromPath` not passing all files with * main.nf is

#!/usr/bin/env nextflow

params.input = ''
tchan = Channel.fromFile(params.input)

tchan.println()

I have a.bam, b.bam and c.bam

I run:

nextflow run main.nf --input *bam

But it only prints a.bam

According to the documentation, wildcard character it is interpreted as a glob path matcher

What am I doing wrong?

Thank you

got it!

I should run:

nextflow run main.nf --input "*bam"

instead of

nextflow run main.nf --input *bam


How to achieve blast results according to the intuitive interpretation of `-max_target_seqs`? Very recently a BLAST parameter -max_target_seqs n got a lot of attention. Instead of the intuitive interpretation (return the best n sequences) the parameters asks blast to return the first n sequences that pass the e-value threshold. This is affecting thousands of workflows and analyses that assumed the intuitive interpretation.

I also found that there is another parameter called -num_alignments. What is the difference? Related question on BioStars seems that is misinterpreting the max_target_seqs parameter. Is num_alignments also misunderstood?

Now the main question. How to actually run blast to give me the best n hits in the database?

Till now I thought that this example will do the job:

blastp -query $QUERY -db $PROTEINS -out $BLASTOUT -evalue 1e-10 -outfmt 6 -num_alignments 5

-- edit --

Now I understood that i can not achieve the desired output only by parameters of blast, but I still don't see an answer to how can I get the best n hits. So I wanted to check out how to get all the results and that's not intuitive either. The parameter -max_target_seqs is set by default to 500. Does it mean that if there are more than 500 significant hits, I don't have a guarantee for the best one? Do I have to specify max_target_seqs to some crazy high number to be sure that I got them all?

The -max_target_seqs is also related to -evalue which does not work as one would think. Looking at the blast news one finds that since version BLAST 2.2.27+ released, one should use the -max_target_seqs:

  4.) The output for reports without separate descriptions and alignments sections (all –outfmt greater than 4) should now use –max_target_seqs to control the output rather than –num_descriptions and –num_alignments.

So it seems that num_alignments is also misunderstood. I came to the conclusion that one cannot expect to get the best n alignments, unless you get all and filter yourself the best ones.


How to filter Ensembl cDNA and ncRNA FASTA files by primary assembly? I'm currently performing differential expression analysis using alignment-free quantification using Kallisto. To do this, I need to create a Kallisto index using Ensembl's cDNA and ncRNA annotations available at the following two links:

ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/cdna/Mus_musculus.GRCm38.cdna.all.fa.gz ftp://ftp.ensembl.org/pub/current_fasta/mus_musculus/ncrna/Mus_musculus.GRCm38.ncrna.fa.gz

However, later on in my analysis I realized that I had multiple Ensembl gene IDs pointing to the same gene. For example, I was gettign results for both ENSMUSG00000026012 and ENSMUSG00000102412 which both map to Cd28. I realized that the reason for this error was due to the fact that the annotations I was using contained haplotypes.

I would like to remove these haplotypes, and to do this I believe that I need to rerun my analysis using a primary assembly. However, Ensembl doesn't have annotated cDNA and ncRNA FASTA files of the primary assembly. How can I filter these files to contain only genes present in the primary assembly? Is there a better way for me to solve this issue?

The simplest method is to download the GTF file for GRCm38 and filter that. You can then use one of the many tools out there (bedtools getfasta (make sure to enforce strandedness and possibly use the --split option, if you happen to use a BED12 file), gtf_to_fasta from tophat2 (it's the only reason to keep tophat2 installed), etc.) to make a transcriptomic fasta file.

Alternatively, Ensembl BioMart can be used to directly download a transcriptomic fasta file using a filter to select only those transcripts on pre-defined chromosomes. Other filters can be applied as well. If you don't already have a GTF file downloaded and one of the previously mentioned packages installed, then this may be the simplest way to get what you want.