Parallel in PLINK for linear association for SNP effects I am doing linear association calculations for SNP effects, I see it take a lot of time to do, does anyone have any experience in parallelizing this job. I have used --threads <num of cores> but it says it only uses 1 thread

here is my code for the calculation in plink:

plink --bfile ukbb_dis --pheno phenotype_dis.pheno --mpheno 4 --allow-no-sex --covar phenotype_dis.pheno --covar-number 5-21 --missing-phenotype -10000000 --linear --adjust --out assoc_SNP_height_linear --threads 4

Here is the current output:

257682 MB RAM detected; reserving 128841 MB for main workspace.
1133273 variants loaded from .bim file.
284516 people (0 males, 0 females, 284516 ambiguous) loaded from .fam.
Ambiguous sex IDs written to assoc_SNP_height_linear.nosex .
284516 phenotype values present after --pheno.
Using 1 thread.
Warning: This run includes BLAS/LAPACK linear algebra operations which
currently disregard the --threads limit.  If this is problematic, you may want
to recompile against single-threaded BLAS/LAPACK.
--covar: 17 out of 21 covariates loaded.
Before main variant filters, 284516 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.992737.
1133273 variants and 284516 people pass filters and QC.
Phenotype data is quantitative.
Writing linear model association results to
assoc_SNP_height_linear.assoc.linear ... 0%

Your help is really appreciated!

Try the latest plink 2.0 alpha build for this.  It’s an incomplete program, but its implementation of —linear is far better than v1.9’s.


What database do gene IDs starting with sp_v2 refer to? I am a student trying to analyze GEO2r datas for one of my courses.
The IDs given in the output are different for different series. I need to convert all of them to a similar format.

In this process I encountered the following type of ID which I don't know its origin:

”ID"    "logFC"
"SP_v2 4634"    "-0.9897758"
"SP_v2 3382"    "-0.8391782"
"SP_v2 4210"    "-1.1693583"
"SP_v2 2117"    "-1.0504727"
"SP_v2 3488"    "-0.9756444"
"SP_v2 1128"    "-0.8289103"
"SP_v2 2735"    "-0.8629999"
...

Each one of the rows represent a single gene.
My question is that what is this ID? The GEO accession is GSE97750.

This is possibly a wild goose chase, but a lot of searching led me to this sample:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM2536852

Which appears to exhibit accessions of the right format. This in turn leads to the platform:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL22166

Which is a custom spotted cDNA array, which the metadata seems to suggest is a human platform.

The full annotation table gives gene symbols and Entrez IDs for most of the probes on the array:

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?view=data&acc=GPL22166&id=53508&db=GeoDb_blob144

To take the example of the first line of your results:

ID          ORF     Entrez gene SEQUENCE
SP_v2 4634  ETNK1   55500       AAAGCAGCTTCATCTTTCAAAATTGATTTGCTCTGGTTTT

The Entrez Gene record for that ID matches up:

https://www.ncbi.nlm.nih.gov/gene/?term=55500


Converter between PDB or mmCIF and MMTF I'd like to test MMTF, a new format for storing biomolecular structures which is promoted by RCSB as a more compact alternative to mmCIF and PDB.

From MMTF FAQ:

  How do I convert a PDBx/mmCIF file to an MMTF file?

  The BioJava library contains methods to read and write PDBx/mmCIF files and MMTF files.

Can I do such a conversion, ideally from command line, but without writing my own Java program?

I wrote a very quick and dirty script to handle conversion between file types using BioJava.

https://github.com/eedlund/Utils/tree/master/BioUtils

Download the jar file here

To run: java -jar BioUtils.jar $FILE $TYPE

where \$FILE is a PDB or mmCIF file you'd like to convert and \$TYPE is the format of the output file [PDB, CIF, MMTF].


Human SNP databases other than dbSNP I'm working on a project to identify variants as being either germline or somatic based on nearby SNPs, and I was wondering if there were any databases of SNPs for Homo sapiens other than dbSNP? I'm not looking for databases that contain additional information (pathogenicity etc.) on SNPs that are already in dbSNP, instead I want to know if there are any other databases that contain SNPs that aren't found in dbSNP.

I've been searching online and so far the only databases I can find are either deprecated or integrated into dbSNP, so any suggestions would be great.

I'm pretty sure there are no public databases of variants that aren't included in dbSNP. That is the variant database so that's the one you want to use.

While there are various other databases around, they tend to use dbSNP rsIDs as well so any variants in those will also be in dbSNP. So while dbSNP is nowhere near complete, of course, it's still the most complete one you're likely to find.


Are the conclusions in "The proximal origin of SARS-CoV-2" legit? I don't have any background in genetics and bioinformatics, so I ask you if you think that the arguments provided in the article The proximal origin of SARS-CoV-2 by Andersen et al. are convincing. In particular:

  While the analyses above suggest that SARS-CoV-2 may bind human ACE2 with high affinity, computational analyses predict that the interaction is not ideal and that the RBD sequence is different from those shown in SARS-CoV to be optimal for receptor binding. Thus, the high-affinity binding of the SARS-CoV-2 spike protein to human ACE2 is most likely the result of natural selection...

I interpret this passage as "if the virus had been engineered, they would have done a better job". Is this interpretation correct? And if it is, does it sound convincing to you as a proof that the virus developed in nature?

The second motivation given in the article is the following:

  Furthermore, if genetic manipulation had been performed, one of the several reverse-genetic systems available for betacoronaviruses would probably have been used

Again, do you think that this exclude the possibility of a human intervention in the creation of the virus?

Thank you for your opinions.

Edit: I am a fan of Occam razor, I know that the scenario in which the virus originated in the wild is by far more likely than the human engineered scenario. I just want to know if, given our current knowledge in genetics, would have been possible for some high skilled researchers to engineered COVID19 (without implying that they did this with any bad purpose)?

Edit 2: I share the concerns of the author of this post, that is, that the findings of the article are merely opinions and some of arguments are misleading. Anyone having a solid technical background can comment on this? China owns Nature magazine’s ass – Debunking “The proximal origin of SARS-CoV-2” claiming COVID-19 definitely wasn’t from a lab

In summary, the authors are saying the complete opposite of "human intervention".

  While the analyses above suggest that SARS-CoV-2 may bind human ACE2
  with high affinity, computational analyses predict that the
  interaction is not ideal and that the RBD sequence is different from
  those shown in SARS-CoV to be optimal for receptor binding.

The interaction with ACE2 is taking place on the Spike protein (S) in both COVID-19 and SARS-Cov. The S protein is very different between COVID-19 and SARS both in terms of sequence and COVID-19 includes an additional furin cleavage site, meaning structurally it is quite distinct. The "computational analysis predicts" is difficult because there is no crystal structure for COVID-19 S protein therefore modelling ligand binding is essentially guess work, particularly as it comprises an additional (furin) cleaved domain. What the authors appear to conclude is that they cannot map ACE2-Spike binding of COVID-19 using an established SARS S protein structure. I assume they are homology modelling SARS-CoV onto COVID-19 S protein sequence data.  In my personal opinion there is no way to homology model an additional furin protease cleavage site (furin is a host protease) so it would be difficult for "good fit" to occur, but it is just my opinion and the authors may have established methods to overcome this notable limitation, e.g. similar phenomena occur in influenza.

  Thus, the high-affinity binding of the SARS-CoV-2 spike protein to
  human ACE2 is most likely the result of natural selection..

Their conclusion of "natural selection" is not strictly accurate and should be "Darwinian positive selection", but its a small issue. What they are saying is Darwinian adaptation has occurred as opposed to "purifying selection, i.e. conservation. Thus they are definitely excluding human intervention because that is not "natural selection"

  Furthermore, if genetic manipulation had been performed, one of the
  several reverse-genetic systems available for betacoronaviruses would
  probably have been used

They appear to be saying that if a synthetic virus was constructed it would be created using an established reverse genetics system, and all those produced to date look nothing like COVID-19.

The phylogenetics argument is always a stronger argument than either of the above (although my opinion isn't entirely objective).  Phylogeneticists have had a long history in opposing conspiracy theories, e.g. polio vaccination, HIV and the list goes on, so it is much more natural territory. An example of such an argument is here which @terdon helpfully pointed out.

To address the questions of @Hans

Question 1

Purifying selection = any evolutionary change is deleterious, so the virus is less able to transmit between humans. The vast majority of mutations are deleterious in classical Darwinian thinking*
Adaptation = the amino acid change occurs to ensure the virus is better able to transmit between humans. In phylogenetics this is called 'positive selection' which is detected at a nucleotide level, this is however a very stringest test.

*, There is something called nearly neutral theory, but thats just complicated

Question 2

What we are talking about is a virus infectious clone system. For a 30kb virus this would not be trival at all, but I agree for 10kb viruses it much easier. However, the phylogenetics are opposed to this, because SARS-CoV-2 is 5% divergent from a virus (RaTG13) isolated a long time before the current epidemic. You can't just engineer a virus that is 5% divergent from a bat virus across the entire genome and never been seen before. It would be a work of genius the likes we have never seen using bioinformatics as yet unheard of, because what would happen is you would encounter endless deleterious mutations. Even if you got over all those, how do you know what fitness traits you are aiming for?

Serial passage
This level of change is doable by serial passage of the virus, but this leads to attenuation of virus (less dangerous) not increased fitness. The yellow fever virus vaccine was created in this way. Even if a virus had been serially passaged we could tell, because we know the mutational patterns of other serially passaged viruses.


How to modify execution orders for Snakemake? My goal is to process two samples, S1 and S2, using two rules, step1, step2.
The Snakemake file is like:

SAMPLES = ['S1', 'S2']

rule all:
  input: expand("{sample}.done", sample = SAMPLES)

rule step1:
  input: "{sample}"
  output: "{sample}.step1.done"
  shell: "touch {output}"

rule step2:
  input: "{sample}.step1.done"
  output: "{sample}.done"
  shell: "touch {output}"

The current execution order is:

step1 on S1
step1 on S2
step2 on S1
step2 on S2

Is there a way to change the order to:

step1 on S1
step2 on S1
step1 on S2
step2 on S2

Basically, I want to complete S1, then S2.
Is there a way to specify the order of executions in Snakemake?

I think you should use the priority directive to give precedence to downstream rule(s).

Here below I give higher priority to step2 so when the first sample has completed step1 snakemake will run step2 on that sample rather than submitting another sample to step1.

SAMPLES = ['S1', 'S2', 'S3']

wildcard_constraints:
    sample= '|'.join([re.escape(x) for x in SAMPLES])

rule all:
    input: expand("{sample}.done", sample = SAMPLES)

rule step1:
    priority: 1
    output: "{sample}.step1.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

rule step2:
    priority: 10
    input: "{sample}.step1.done"
    output: "{sample}.done"
    shell:
        r"""
        sleep 5
        touch {output}
        """

(Note that without wildcard_constraints your/this code throws AmbiguousRuleException)


How can I get plant Genome size, total number of genes, EST, and chromosomes number How can I get plant Genome size, total number of genes, EST and chromosome number from phytozom and plantDB?

Ensembl plants has a REST API, you might be able to get all these details (except possibly EST numbers) there: http://ensemblgenomes.org/info/access/rest

For example this query gives you the Arabidopsis genome size: http://rest.ensemblgenomes.org/info/genomes/arabidopsis_thaliana?content-type=application/json


How to check if a given gene is expressed in a group of microarray samples if I do not have control group to compare with? I have a microarray gene expression dataset consisting of placenta samples. I want to check whether these placenta samples are mixed with maternal decidua tissue. I have marker genes for decidua samples SCARA5, PIGF, SLPI. How to check whether these genes are expressed in placenta samples? In other words, how to check if a gene is up or down regulated without control group?

You ask about which "genes are expressed" and then you mention "if a gene is up or down regulated". These are different, and given your application I think what you actually want to know if these marker genes are expressed. This is not a question of up/down-regulation relative to a control, so you do not need control samples for this.
It seems most of your datasets are Affymetrix, and there are R/Bioconductor packages that will give you the probability a gene is expressed. One such I found is frma (https://bioconductor.org/packages/release/bioc/html/frma.html), which returns an object "weight" that is "a vector of weights which roughly correspond to the probability of expression for each gene." I haven't worked with microarrays in a while, but I recall many other R/Bioconductor packages offering this functionality.
A detection p-value is calculated for Illumina array probes in Illumina's native BeadStudio software.


taxon exclude list for searching local blast database using blastn I am looking for a solution to exlude certain entries when searching a local blast nt database (with blastn), specifically the sequences from uncultured / environmental samples, ideally using their taxon ids.
This should ideally replicate the results from the BLAST web interface, when checking the option "Exclude Uncultured/environmental sample sequences", see screenshot below.

Three approaches that I am aware of so far are:

(1) grep taxon IDs for any taxa with "uncultured" or "environmental" in their scientific name from names.dmp and fetch the GIs belonging to these taxa using the E-utils. Then use the -negative_gilist parameter for blastn. However, GIs are discontinued for future records, so this method should be replaced by Acc.Ver identifiers.

(2) locally format the blast database after prepocessing the nt file by myself and exclude records belonging to the taxa list as obtained in (1).

(3) Postprocess the blast output and remove the relevant lines, but this approach might remove a lot of alignments, whereas (1) and (2) would always generate the desired number of output alignments.

So, I am looking for other suggestions, preferably without preprocessing the nt file.
Ideally, one would like to have a -negative_taxalist option for blastn, which would resemble what is possible on the NCBI BLAST web interface, where one can exclude taxa, see screenshot below.

Thanks!

Hurray, the solution seems already to be implemented in the new BLAST database format (v5), which is usable with blast+ versions >= 2.8.0.

It is possible to specify an exclusion as well as inclusion list of taxon ids, which, however, have to be on species level.

Species level taxon ids can be obtained from higher level taxon ids using the script get_species_taxids.sh, e.g.: get_species_taxids.sh -n Enterobacterales

See docs here: https://ftp.ncbi.nlm.nih.gov/blast/db/v5/blastdbv5.pdf

Edit:
For creating taxon ID lists locally from the NCBI taxonomy files, the program subtree can be used. See example therein for compiling a list of taxon IDs from environmental or unclassified taxa.


Calculating the mortality rate of a pandemic, e.g. coronavirus 2019-nCov? When calculating how lethal a virus outbreak is, I've noticed that most sources use

CONFIRMED_DEATHS/CONFIRMED_INFECTIONS

However, given that there are only 2 possible outcomes; death or recovery, wouldn't it be more logical to use:

CONFIRMED_DEATHS/(CONFIRMED_DEATHS + CONFIRMED_RECOVERIES)

Which formula is applicable to which situations?

This may not strike most as a bioinformatics, but getting the key clinical outcome is essential in understanding the molecular basis of pathogenicity.

I think the mortality rate is over-reported. This is not to say the situation of 2019-nCov is not serious - it is very serious.

The two essential factors missing in your equation for 2019-nCov are:

Age. Mortality is likely skewed towards the old and very young
Seroprevalence. To identify the number of asymptomatic patients.

Old people have a weakened immune system (e.g. reduced T-cell response). The very young have an underdeveloped immune system. The problem is hospital cases are skewed towards those with severe infection, which is likely the old and very young. In 2019-nCov the first infecteds were middle aged men, but likely to be restricted to those exposed to the "Wuhan meat market".

Problem
Asymptomatic "spreaders" are known in influenza and are suspected in 2019-nCov. Asymptomatics will inflate the mortality rate (because we're not detecting them).

Solution(s)

One solution is to conduct seroprevalence surveys. Serological surveys screening for antibodies (IgG) against 2019-nCov will identify asymptomatic patients, using the first formula you stated. However, would you seriously go around the general population in a 2019-nCov epidemic taking blood samples? Errhhh no.
The bio-statistics approach is to use a catalytic model and critical to this is the age of the patient. Catalytic models are a priori models that incorporates the age-prevelance profile of infection. By introducing the incidence and age of hospitalisation of e.g. the old and the young, the rest of the distribution, i.e. asymptomatics can be calculated by regression. This estimates the total population of infecteds and again you use equation 1 to calculate the mortality rate recorded by the hospital.
Mathematical models based on transmission dynamics. The problem here is the rate of transmission is unknown.

I would assume the hospital is the place virus-induced mortality is reasonably accurate.

The second formula you described works very well for Ebolavirus, no-one has ever heard of asymptomatic Ebolavirus. It would also work well for another betacoronavirus Middle Eastern respiratory syndrome. The key problem with 2019-nCov is we don't know the number of asymptomatic patients or patients with such mild infections they never report illness.

Summary
In summary, I think the mortality is being over-reported if it is based on your equation 2. Equation 1, however, will also over-report the mortality rate.


Get a list of PDB ids from a list of sequence protein number? I am working on a project about phosphorylation, and I am creating a database of some proteins. I have a list of protein accession ids for which I need to get the corresponding PDB file for each one to map it to PSM matrix. I want to get the PDB file using Python or BioPython tools. Any suggestions will be appreciated.

A typical Uniprot URL would look like the example here:

https://www.uniprot.org/uniprot/?query=ID:P26378&format=tab&columns=id%2Cdatabase(PDB)

The response will be:

  Entry Cross-reference (PDB)

  P26378    1FXL;1G2E;

We want our query to include the protein ID, the output format should be tab separated and the only contain the columns id and database (PDB).

To get the PDB IDs you would need to split the response, take the last line, split this line after each tab and split it again after each semi-colon.

Using Python's Requests library the whole code would be.

import requests
import time

uniprot_ids = ['P26378', 'O35433', 'Q02910']
url = 'https://www.uniprot.org/uniprot/'

protein_to_pdb = {}
for protein in uniprot_ids:
    params = {
        'format': 'tab',
        'query': 'ID:{}'.format(protein),
        'columns': 'id,database(PDB)'
    }
    contact = ""  # Please set your email address here.
    headers = {'User-Agent': 'Python {}'.format(contact)}
    r = requests.get(url, params=params, headers=headers)

    protein_to_pdb[protein] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')
    protein_to_pdb[protein].pop(-1)
    time.sleep(1)  # be respectful and don't overwhelm the server with requests

print(protein_to_pdb)

To accelerate things, you could store your protein_to_pdb dict in a file and only retrieve the new ones
Checking if the returned ID matches your input ID is always a good idea


Difference between BWA-backtrack and BWA-MEM Many of my colleagues recommend I use BWA-MEM instead of regular old BWA. The problem is I don't understand why and reading the BWA man page doesn't seem to help the matter.

What is the difference between BWA and BWA-MEM? And, in which instances would you employ one over the other?

TL;DR:

BWA-backtrack is based on backtracking.
This approach is appropriate only when the dissimilarity between the reads and the reference is low,
or when you want to find all best hits or enumerate all possible alignments up to a specified
number of errors.

In all other situations, BWA-MEM is preferable as it can,
thanks to its sophisticated strategy based on maximum exact matches,
deal with errors better and also automatically switch between local and
global alignment modes.

Long description:

I would like to provide some algorithmic insight since I believe that it
can be very useful in this case.
Both BWA-backtrack and BWA-MEM use the same indexing strategy
(heavily relying on BWT-index), but the actual algorithms are quite different.

BWT-index (and also other full-text indexes such as
suffix arrays or suffix trees)
can easily find exact matches (imagine Ctrl+F-like search in a text editor),
but any differences between the read and the reference,
such as sequencing errors or genomic variants,
make the situation complicated. One then needs to
somehow transform inexact matching to exact matching,
and all three BWA mappers (note that there exists also BWA-SW, but it is deprecated)
use quite different strategies.

BWA-backtrack looks for substrings of the reference, which would be similar
to the entire read (end-to-end) using an algorithm
called backtracking.
First, it searches occurences of the read without any "corrections".
If nothing found, it consideres all possible single edits; then two edits, etc.
To make mapping efficient, one usually wants to stop with the first
found alignment as this would be the best one.
When required, it is also possible to find the other equally good alignments
or to enumerate all alignments up to some edit distance or withing some divergance rate
(see the -N option of BWA-backtrack).

It turns out that
the time required for finding an alignment
can be exponential in the number of errors, which is probably the main problem of
backtracking-based approaches.
To prevent huge overheads due to dissimilar reads,
one needs to limit the number of allowed errors to some reasonable number
(see maxDiff in the BWA man page)
and consider the other reads unaligned.
In the case of BWA-backtrack,
the minimum required identity level is ~97% with the default options (see the -n option).

In fact, the algorithm is more complicated and uses various heuristics such as
seed-and-extend or Z-dropoff in order to make the computation fast enough (at the price of lower accuracy). If you are interested
in more details, all these tricks are well described in the paper.

BWA-MEM uses quite a different strategy. It detects long exact matches
between the read and the reference, and then chains them into local or global alignments,
based on what is more appropriate in that specific case.
Such an automatic local-global switching can be very powerful and
BWA-MEM works well with various types of data
(short reads, long reads, low error rates, high error rates, etc.).


Error correction within the long read I am going to get some data from plasmid sequencing to identify SNPs on the plasmids. What is done in the lab is the following:

The plasmids are purified by size.
We amplify the plasmids using the phi29 polymerase. The polymerase will go through the plasmid multiple times. Hence we get the same sequence concatenated multiple times. My question is related to this step:
Finally, We sequence it using Oxford Nanopore.

My question is:

On step two, I wrote that "we get the same sequence concatenated multiple times". For me, this a potential source of information to correct the base calls prior to aligning them to the reference. Since you have the same sequence multiple times (concatemers).

What I would like to know is if there is a tool that uses the information of the concatemers to improve the base calls. I have tried to find some methods but found none. I know that PacBio has a similar flavour using the "circular consensus calling", but I have not found any methodological explanation.

thanks!

Yes, there is a tool to do this called R2C2 by the Vollmers lab at UC Santa Cruz.

https://www.biorxiv.org/content/early/2018/06/04/338020


Assessing PyMol sequence alignment object I use cealign in PyMol for structure alignment. Instead of visualization, I want to return the alignment object to my python script for further analysis. Is there a function to return the object without I/O to a clustalw file?
edit:
By alignment object, I meant the sequence alignment in cealign.

In the application, unfortunately the only option is doing

cealign ,, object =
Menu > Export Alignment... on your cealign  object
and converting the Clustal format using one of the millions of tools online

get_fastastr
As you know cealign adds gaps, but cmd.get_fastastr('') strips them. There is no option to keep the gaps —annoyingly especially with missing loop info!
The help has this to say:
 get_fastastr(selection='all', state=-1, quiet=1, key='', _self=<module 'pymol.cmd' from '/Users/matteo/miniconda3/lib/python3.7/site-packages/pymol/cmd.py'>)

Correction and dissection

So I had a quick check to see if there is an alternative route to get the sequence string with gaps. These are my observations using the Python3 module, sorry if they aren't polished up.

In the comments I said something incorrect. cmd.cealign returns a dictionary with the 4x4 rototranslational matrix —no need for an object. Using Python3 and two TIM barrels as an example:
import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1TIM', 'tim')
    pymol.cmd.remove('not chain A')
    pymol.cmd.fetch('5AHE', 'hisA')
    x = pymol.cmd.cealign('tim', 'hisA')
    print(x)
# {'alignment_length': 184, 'RMSD': 3.8800425577361755, 'rotation_matrix': [...]}

While the object does not contain any useful info from PyMol directly*:
    ...
    pymol.cmd.cealign('tim', 'hisA', object ='rototrans')
    pymol.cmd.get_session('rototrans')['names'][0][5]
    # useful info is generally in ['names'][0][5]

(*) The code is actually wrapped C++ code, so the object rototrans is in the C++ layer and can only be access by session or possibly by some functions in cmd.exporting. In fact, it can be exported to an aln string (ClustalW).
    aln = pymol.cmd.exporting.get_alnstr('rototrans')

In the PyMOL application, doing the cealign with an object (cmd.cealign('tim', 'hisA', object ='rototrans'), allows the  Menu>Export Alignment...>'rototrans' which is saved as a ClustalW file, as you mention. Namely:
CLUSTAL

tim          APRKFFVGGN----WK-------MNGKRKSLGELIHTLDGAKLSADTEVVCGAPSI------Y
hisA         MIIPALDLIDGTVVRLHQGDYARQRDYGN---DPLPRLQDYAAQGAGVLHLVDLTGAKDPAKR
                           ..       ......   .....*..................      .

tim          LDFARQK---LDAKIGVAAQNCYKVPKGAFTGEISPAMIKDI-GAAWVILGHSERRHVFGESD
hisA         QIPLIKTLVAGVNVPVQVGGGVRTEE--------DVAALLKAGVARVVIGS---------TAV
             .......   ................        ..*..... .*..**..         ...

tim          ELIGQKVAHALAEG---LGVIACIGEKLDEREAGIT-----------EKVVFQETKAIADNVK
hisA         KSPDV-VKGWFERFGAQALVLALDV---RIDEHGTKQVAVSGWQENSGVSLEQLVETYLPVGL
             ..... *.......   ..*.*...   ...*.*..           .....*..........

tim          DWSKVVLAYEPVWAIGTGKTATP------QQAQEVHEKLRGWLKTHVSDAVAVQSRIIYGGSV
hisA         ----KHVLCTDI----------SRDGTLAGSNVSLYEEVCARYP---------QIAFQSSGGI
                 ........          .      .......*.......         *......*..

tim          TGGNCKELASQHDVDGFLVGGASLKPEFVDIINAKH------------
hisA         GDIDDIAALRGTGVRGVIVGRALLEGKFTVKEAIQ?WQNVKGHHHHHH
             .............*.*..**.*.*...*..

However, there is no method in cmd.exporting that gives out just the sequence with gaps.
Whereas the cmd.exporting has the following... none applicable.
 print(dir(pymol.cmd.exporting))
['DEFAULT_ERROR', 'DEFAULT_SUCCESS', 'QuietException', 'SDF', 'SDFRec', 'Shortcut', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_cmd', '_feedback', '_get_dump_str', '_get_mtl_obj', '_raising', '_resn_to_aa', '_session_convert_legacy', '_unit2px', 'absolute_import', 'assign_atom_types', 'cPickle', 'cache', 'cache_action_dict', 'cache_action_sc', 'cmd', 'copy', 'copy_image', 'export_coords', 'fb_mask', 'fb_module', 'get_alnstr', 'get_bytes', 'get_ccp4str', 'get_cifstr', 'get_fastastr', 'get_maestr', 'get_mol2str', 'get_pdbstr', 'get_pqrstr', 'get_psestr', 'get_sdfstr', 'get_session', 'get_str', 'get_xyzstr', 'io', 'is_dict', 'is_error', 'is_list', 'is_ok', 'is_tuple', 'loadable', 'lock', 'multifilenamegen', 'multifilesave', 'multisave', 'os', 'png', 'print_function', 'pymol', 're', 'save', 'savefunctions', 'selector', 'sys', 'thread', 'traceback', 'unlock']


VCF-merge fails due to tabix not producing .tbi files I'm trying to use vcf-merge to combine 2 exome capture vcf files (~250K and ~330K in size) before trying it on all 96 samples. I'd appreciate any advice on the best way to do that! I've detailed what I've tried below. My issue seems to be with using tabix to convert the files to .tbi format.

Step 1: BGZIP

So far, I've zipped the files without issue:

bgzip sample1.vcf
bgzip sample2.vcf

which produces:

sample1.vcf.gz
sample2.vcf.gz

Step 2: TABIX

When I try to use this command:tabix -h -p vcf sample1.vcf.gz the stderr is: Region 536999277..536999278 cannot be stored in a tbi index. Try using a csi index with min_shift = 14, n_lvls >= 6. tbx_index_build failed: sample1.vcf.gz

Using the -C option which works:

tabix -C -h -p vcf sample1.vcf.gz
tabix -C -h -p vcf sample2.vcf.gz

which produces:
    sample1.vcf.gz.csi
    sample2.vcf.gz.csi

Step 3: VCF-MERGE

When I use this command: vcf-merge sample1.vcf.gz.csi sample2.vcf.gz.csi > out.vcf.gz the merge fails, and I get this stderr:

Broken VCF header, no column names? at /usr/share/perl5/Vcf.pm line 172, <__ANONIO__> line 1. Vcf::throw(Vcf4_2=HASH(0x5645a760bc38), "Broken VCF header, no column names?") called at /usr/share/perl5/Vcf.pm line 867 VcfReader::_read_column_names(Vcf4_2=HASH(0x5645a760bc38)) called at /usr/share/perl5/Vcf.pm line 602 VcfReader::parse_header(Vcf4_2=HASH(0x5645a760bc38)) called at /usr/bin/vcf-merge line 183 main::init_cols(HASH(0x5645a761f438), Vcf4_2=HASH(0x5645a760b248)) called at /usr/bin/vcf-merge line 279 main::merge_vcf_files(HASH(0x5645a761f438)) called at /usr/bin/vcf-merge line 12

If I exclude the tabix files, the merge still fails and the stderr says: The column names not tab-separated? Could not load .tbi index of sample1.vcf.gz. The command exited with an error. Is the file tabix indexed?

Any suggestions on how I can correctly merge these files?

You should go with the CSI index as it seems like the BAI indexes cannot handle chromosomes longer than 2^29.

Your problem is that you're trying to merge the index files using vcf-merge. You should merge the VCF files, not the .csi index files.


programmatic secondary structure prediction for >36-mer DNA oligonucleotides I'm writing a tool to automate Sanger sequencing primer design for a production lab that uses a universal-tail chemistry Sanger sequencing to verify NGS results. Essentially, the template DNA is amplification using primers with universal 5' tails (M13F/R), and then the "BigDye" reaction is performed with an M13F and M13R primers in separate wells.

I'm able to design logical primers for each region in a request using Primer3 (v2.5.0), but then I'd like to calculate folding thermodynamics for the M13-appended primers to determine the best primer partner for the 5' M13 appendages.

In other words, in the last step, I would like to know if I should use:

  M13F-forward_primer + M13R-reverse_primer

or

  M13R-forward_primer + M13F-reverse_primer

The problem is that Primer3 won't fold an oligo > 36 nt, since the algorithm for Tm is diminishingly accurate for longer oligos. However, I believe thermodynamic predictions are valuable, even if only relatively, i.e. comparing one 5' appendage to another 5' appendage on the same allele-specific sequence, so I'd like the most painless answer, even if it isn't accurate in absolute terms.

I'm aware that UNAFold will do this, but it seems like extreme overkill for what I need to do.

For lack of a better solution, I've implemented a workaround wherein I pair each allele-specific primer with either M13F or M13R and submit the pair as a check_primer task. In other words, the left primer is one of the allele-specific primers that primer3 designed, and the right primer is either M13F, or M13R. I then choose the M13 appendage orientation based on the PRIMER_PAIR_0_COMPL_ANY_TH output values for the two possible orientations.

I would like to know: Is there a better solution? Primer3 docs say that the Tm calculation (and not necessarily folding thermodynamic calculations) is inaccurate for 37+ mers. Is there a way to force Primer3 to fold a sequence > 36 nt? If not, is there a simple tool for folding >36-mers?

Follow-Up: The workaround described in my OP was inadequate and UNAFold would cost $9k (which would take months for my employer to approve - don't get me started), so I tried @MichaelG's solution below. It turned out to work quite well. I'm computing the maximum Tm value for all fragments from every oligo, and that number agrees quite well with the browser-based UNAFold output.

If you really insisted on using Primer3 you could perform a 36bp sliding window analysis. You would submit your primers in 36bp 'chunks' with a step size of say 3 bp. You would need to "circularise" your primer within the sliding window to assess the "beginning" and the "end" and remove "artifactual" hairpins. Its doable via a script and you would automate the web upload and data extraction via a "digger/spider" or else just use a download equivalent (I've used a standalone). There would be a "blindspot" using this approach, which you could overcome if you wanted.

So summary, its complicated and I think your suggestion of using UNAFold is a better idea.


Seurat installation fails When I am trying to install Seurat in R markdown file:

install.packages("Seurat")

it fails with the following errors:

Installing package into ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3’
(as ‘lib’ is unspecified)
Warning in install.packages :
  dependency ‘proxy’ is not available
also installing the dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’

trying URL 'https://cran.rstudio.com/src/contrib/igraph_1.2.1.tar.gz'

Content type 'unknown' length 2703426 bytes (2.6 MB)

downloaded 2.6 MB

trying URL 'https://cran.rstudio.com/src/contrib/dtw_1.18-1.tar.gz'

Content type 'unknown' length 849256 bytes (829 KB) downloaded 829 KB trying URL 'https://cran.rstudio.com/src/contrib/diffusionMap_1.1-0.tar.gz'

Content type 'unknown' length 39130 bytes (38 KB)
downloaded 38 KB

trying URL 'https://cran.rstudio.com/src/contrib/Seurat_2.3.0.tar.gz'

Content type 'unknown' length 846065 bytes (826 KB)

downloaded 826 KB
* installing *source* package ‘igraph’ ...
** package ‘igraph’ successfully unpacked and MD5 sums checked
installing to /sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs
** R
** demo
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded

Error in dyn.load(file, DLLpath = DLLpath, ...) :    unable to load shared object '/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs/igraph.so': libicui18n.so.58: cannot open shared object file: No such file or directory

Error: loading failed Execution halted ERROR: loading failed
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph’ Warning in install.packages :   installation of package ‘igraph’ had non-zero exit status

ERROR: dependency ‘proxy’ is not available for package ‘dtw’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/dtw’ Warning in install.packages :   installation of package ‘dtw’ had non-zero exit status

ERROR: dependency ‘igraph’ is not available for package ‘diffusionMap’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/diffusionMap’ Warning in install.packages :   installation of package ‘diffusionMap’ had non-zero exit status

ERROR: dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’ are not available for package ‘Seurat’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/Seurat’ Warning in install.packages :   installation of package ‘Seurat’ had non-zero exit status

How could it be fixed? What to try?

Operating system is CentOS Linux release 7.3.1611 (Core)

  Update

I also tried development version of Seurat but it failed too:
http://satijalab.org/seurat/devinstall.html

  Update

I tried installing Seurat using conda install -c bioconda r-seurat and it is not working. It is giving the error:

  -bash-4.2$ conda install -c bioconda r-seurat
  Fetching package metadata .............
  Solving package specifications:

  PackageNotFoundError: Packages missing in current channels:

r-seurat -> r-cowplot
r-seurat -> r-diffusionmap
r-seurat -> r-dtw
r-seurat -> r-ggridges
r-seurat -> r-ica
r-seurat -> r-pbapply
r-seurat -> r-ranger
r-seurat -> r-rcppprogress
r-seurat -> r-sdmtools
r-seurat -> r-tclust

  We have searched for the packages in the following channels:

https://conda.anaconda.org/bioconda/linux-64
https://conda.anaconda.org/bioconda/noarch
https://repo.continuum.io/pkgs/main/linux-64
https://repo.continuum.io/pkgs/main/noarch
https://repo.continuum.io/pkgs/free/linux-64
https://repo.continuum.io/pkgs/free/noarch
https://repo.continuum.io/pkgs/r/linux-64
https://repo.continuum.io/pkgs/r/noarch
https://repo.continuum.io/pkgs/pro/linux-64
https://repo.continuum.io/pkgs/pro/noarch

As @terdon mentioned in the comments, you need to sudo yum install libicu on the relevant nodes (or have your admin do that for you). Alternatively, seurat can be installed via conda, which means you don't need root access.


Strange p-value histogram for differential gene expression analysis I'm trying to perform pretty standard differential expression analysis using RNA-seq data.

I've used Kallisto to perform RNA quantification and am using Sleuth to perform the differential expression analysis.

I have 2 conditions that I'd like to compare, which is stored in the metadata:

> metadata
    sample cell replicate batch                   path
1:  Base-1 Base         1   1   /path/to/Batch1/Base-1
2:  Base-2 Base         2   2   /path/to/Batch2/Base-2
3:  Base-3 Base         3   2   /path/to/Batch2/Base-3
4:  Test-1 Test         1   2   /path/to/Batch2/Test-1
5:  Test-2 Test         2   2   /path/to/Batch2/Test-2
6:  Test-3 Test         3   3   /path/to/Batch3/Test-3

I'm attempting to test for differentially expressed genes between the Base and Test cells while controlling for batch (the sequencing batch). The R code to do this is adapted from Sleuth's manual:

sleuth_obj = sleuth_prep(
    metadata,
    extra_bootstrap_summary = TRUE,
    num_cores = 1,
    target_mapping = transcripts_genes, #ENSEMBL information mapping transcripts to genes
    aggregation_column = "ens_gene"
)

# smooth raw Kallisto abundances
sleuth_obj = sleuth_fit(sleuth_obj, ~cell + batch, "full")

# fitting error measurements on reduced model
sleuth_obj = sleuth_fit(sleuth_obj, ~batch, "reduced")

# compare models to calculate differentially expressed transcripts
sleuth_obj = sleuth_lrt(sleuth_obj, "reduced", "full")

# extract results
sleuth_table = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE
))
sleuth_table_tx = as.data.table(sleuth_results(
    sleuth_obj,
    'reduced:full',
    'lrt',
    show_all = FALSE,
    pval_aggregate = FALSE
))

To check the quality of the null hypothesis, I've checked the p-value histograms of the gene-level and transcript-level analyses, and I get the following:

Neither of these plots show well-behaved p-values, which makes me question the validity of the reduced model. How can I go about adjusting my analysis to get more reliable results?

Edit following from @swbarnes2's answer:

I've used fdrtool to create an empirical null model of Sleuth's likelihood ratio test.

library(fdrtool)
local_fdr_txs = fdrtool(so$tests$lrt$`reduced:full`$test_stat, statistic="normal")

Which gives the following plots:

Plotting the distribution of these p-values produces:

This still has a hump, but looks more anti-conservative than before.

My remaining question is about whether I can use fdrtool to model the null hypothesis as normal, since I'm not certain about the distribution for Sleuth's likelihood ratio test

Try the protocol here for generating a different (smaller) dispersion estimate

https://www.huber.embl.de/users/klaus/Teaching/DESeq2Predoc2014.html#inspection-and-correction-of-pvalues


removing indels in plink I have data in Plink format. I removed Indels using --snps-only command in Plink, but still some snps are as below:

rs58119544;rs5746975    0   16915280    T   A

I wonder, how i can get rid of these snps? the correct should be single snps without this (;) sign. Any idea to fix it please?

Use grep to get SNPs that have ";" in the name:

>$ grep ";" myMap.txt | cut -f1 -d" " > myBadSNPs.txt

Then use plink with --exclude :

plink --file myFile
--exclude myBadSNPs.txt
--out myFileFiltered


Sophisticated algorithms in biological systems I am looking for examples of (known) algorithms/data structures in biological systems. My current understanding is that these don't exist - or haven't been found - beyond a relatively low level of complexity, which would seem like strong evidence to me that evolution isn't able to come up with something as complicated (and multi-step) as, e.g., the Schönhage-Strassen algorithm.

I heard of Bloom filters for recognizing smells in fruit flies (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6304992/).

There are many examples of biological systems doing computation, according to my understanding:

Regulatory networks encode decision trees of significant complexity: wikipedia.
Alan Turing did significant work on organismal morphogenesis showing how simple chemical patterns of inducers can lead to complex outcomes, "Turing patterns".
All of organismal development, cell division, etc. are tightly regulated procedures of gene regulation and morphogenesis.
The inspiration for neural networks is brains.


How can longest isoforms (per gene) be extracted from a FASTA file? Is there a convenient way to extract the longest isoforms from a transcriptome fasta file? I had found some scripts on biostars but none are functional and I'm having difficulty getting them to work.

I'm aware that the longest isoforms aren't necessarily 'the best' but it will suit my purposes.

The fasta was generated via Augustus. Here is what the fasta file looks like currently (sequence shortened to save space)

>Doug_NoIndex_L005_R1_001_contig_2.g7.t1
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca
ttaatataaccaaaaacctagaaattctagccttcgatgatgttgcagctgcggttcttgaagaagaaagtcggcgcaagaacaaagaagatagaccg
>Doug_NoIndex_L005_R1_001_contig_2.g7.t2
atggggcataacatagagactggtgaacgtgctgaaattctacttcaaagtctacctgattcgtatgatcaactcatca

The format is as such:

Gene 1 isoform 1
Gene 1 isoform 2
Gene 2 isoform 1
Gene 2 isoform 2

and so forth. There are several genes that have more than one pair of isoforms (up to 3 or 4). There are roughly 80,000 total transcripts, probably 25,000 genes. I would like to extract the single longest isoform for each gene.

While the solution from https://bioinformatics.stackexchange.com/users/96/daniel-standage should work (after adjusting for possible python3 incompatibility), the following is a shorter and less memory demanding method that uses biopython:

#!/usr/bin/env python
from Bio import SeqIO
import sys

lastGene = None
longest = (None, None)
for rec in SeqIO.parse(sys.argv[1], "fasta"):
    gene = ".".join(rec.id.split(".")[:-1])
    l = len(rec)
    if lastGene is not None:
        if gene == lastGene:
            if longest[0] < l:
                longest = (l, rec)
        else:
            lastGene = gene
            SeqIO.write(longest[1], sys.stdout, "fasta")
            longest = (l, rec)
    else:
        lastGene = gene
        longest = (l, rec)
SeqIO.write(longest[1], sys.stdout, "fasta")

If you saved this as filter.py, then filter.py original.fa > subset.fa would be the command to use.


Strange per sequence GC content results I would be grateful if someone could take a quick look at these FASTQC results. This is rna-seq paired-end data. From the FASTQC manual, an unusual distribution seems to be suggestive of contamination and a shift in the curve is suggestive of a systematic bias. GC content distribution both pre-alignment and post-alignment are strange.

Samples are paired end, strand specific and % of mapped reads is above 95% for all the samples. There is no adapter content also.

What could be the problem?

Although generally I recommend people to be wary about the warnings that FastQC gives (it tends to be overly paranoid), the GC content graphs here do look odd. It's good that you've had a look at it with another program to confirm the observation. They should have a distribution that is fairly close to normal, a bit like this:

Without more information about your sample, I can't tell whether or not this is expected from your input. This looks like fairly substantial contamination, but could also be caused by biased selection during sample preparation (e.g. low quality sample, low yield, PCR amplification bias), or just looking at a multi-organism sample rather than a single-organism sample.

As a contamination check, I'd recommend running your RNASeq reads through a kmer-based classification tool (e.g. OneCodex, Centrifuge, BBSketch, or Sourmash). This should tell you the most likely place that contaminant sequences came from.


Plink 1.9 --merge-list changing order of fam file I am doing some whole-genome sequence analysis and to speed things up split the data into chromosomes and also into regions within those chromosomes. (So I have my dataset split into ~300ish files) I am doing some IBD analysis with plink files and to do this IBD analysis I am combining the files into one big plink file.

I ran into issues today where when using plink 1.9 and the --bfile and --merge-list commands, the output fam file is different from the input fam file.  My dataset has related individuals so the fam file is very important.  Also its is important to note that when using binary files the order in the fam file needs to be exactly the same as the order the individuals appear in the other files (I learned this the hard way haha).

Does anyone know why plink 1.9 changes the order in the fam file when merging plink files, and how to fix this?  It works with plink 1.07 but with 1.9 I get a different fam file.

Performing "diff" on input and output fam files shows they are different:

plink=$mydir/bin/plink_1_9/plink

$plink --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

Using the same code but using plink-1.07 creates matching fam files as expected"

plink1=bin/plink-1.07-x86_64/plink
$plink1 --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

The merge_list.txt file is in this format:

$indir/chr22/286/$prefix.chr.22.region.286.qc.bed $indir/chr22/286/$prefix.chr.22.region.286.qc.bim $indir/chr22/286/$prefix.chr.22.region.286.qc.fam
$indir/chr22/287/$prefix.chr.22.region.287.qc.bed $indir/chr22/287/$prefix.chr.22.region.287.qc.bim $indir/chr22/287/$prefix.chr.22.region.287.qc.fam
$indir/chr22/288/$prefix.chr.22.region.288.qc.bed $indir/chr22/288/$prefix.chr.22.region.288.qc.bim $indir/chr22/288/$prefix.chr.22.region.288.qc.fam

Can someone explain to me why this could be happening and if there is a way to fix it? Thanks so much!!

(Note, I know I can manually replace the fam file with the correct one with cp or mv, but I would like to know how to handle this with plink 1.9, if possible, as I am afraid it could be indicative of an underlying error in the merge process :) )

Thanks Christopher Chang for the great answer via the plink2 google group! (See here)

His answer was as follows:

"Yes, plink 1.9 did change the default merged sample order.  However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge.  (You can also use --indiv-sort + --make-bed at any time to change the sample order to whatever you want.)

As for why the default was changed, take a look at the new default order: it will probably make a lot of sense.  It’s “natural-sorted” which tends to make it easy for humans to manually search through regardless of what mix of letters and numbers are used in the IDs."

See also the Plink 1.9 documentation for --indiv-sort which makes the point that this behavior is default only for merges--in other commands Plink 1.9 uses the same ordering as in Plink 1.07, which is the order the samples were in.


Parsing MSigDB Supplementary Collections MSigDB recently added Supplementary Collections. The formatting is similar to the original MSigDB gene sets, but there is an additional metadata.txt file (link).

It looks like (condensed for easier reading):

STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...
STANDARD_NAME   Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron
ORGANISM    Homo sapiens
ORGAN_SYSTEM    Central Nervous System
PMID    29867213
PUBLICATION_TITLE   Spatial transcriptomic survey of human embryonic cerebral cortex by single-cell RNA-seq analysis
AUTHORS Fan X,Dong J,Zhong S,Wei Y,Wu Q,Yan L,Yong J,Sun L,Wang X,Zhao Y,Wang W,Yan J,Wang X,Qiao J,Tang F,
GEOID   GSE103723
...

I could not find any description of the formatting. Is this a standard format? Is there a parser for it (preferably in R)? It looks like I can just transpose each record to get it in a table format, but I don't know if I can trust all records to be consistent.

If you want to separate each set, you can do something like this

dat <- readLines("~/Desktop/scsig.v1.0.metadata.txt")
dat <- gsub("\t"," ",dat)

# get the indeces of each type
standardName <- grep("STANDARD_NAME",dat)
organism <- grep("ORGANISM",dat)
organSystem <- grep("ORGAN_SYSTEM",dat)
pmid <- grep("PMID",dat)

# return the values
standardName <- gsub(pattern = "STANDARD_NAME",replacement = "",dat[standardName])
organism <- gsub(pattern = "ORGANISM",replacement = "",dat[organism])
organSystem <- gsub(pattern = "ORGAN_SYSTEM",replacement = "",dat[organSystem])
pmid <- gsub(pattern = "PMID",replacement = "",dat[pmid])

df <- data.frame(standardName,organism,organSystem,pmid,stringsAsFactors = FALSE)

 head(df)
                                     standardName      organism             organSystem      pmid
1      Fan_Embryonic_CTX_Big_Groups_Cajal_Retzius  Homo sapiens  Central Nervous System  29867213
2  Fan_Embryonic_CTX_Big_Groups_Brain_Endothelial  Homo sapiens  Central Nervous System  29867213
3  Fan_Embryonic_CTX_Big_Groups_Excitatory_Neuron  Homo sapiens  Central Nervous System  29867213
4              Fan_Embryonic_CTX_Big_Groups_Glial  Homo sapiens  Central Nervous System  29867213
5       Fan_Embryonic_CTX_Big_Groups_Brain_Immune  Homo sapiens  Central Nervous System  29867213
6          Fan_Embryonic_CTX_Big_Groups_Microglia  Homo sapiens  Central Nervous System  29867213

This method is very verbose, but it's a start and it can be easily improved with a loop.

 #---- EDIT

Not knowing how the format works, there may be two issues:

1) gene sets with missing values only

STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
ORGANISM
ORGAN_SYSTEM
PMID 31292543

2) gene sets with missing value and type

STANDARD_NAME Aizarani_Liver_C33_Stellate_cells_2
PMID 31292543

Case 1 is not an issue, as the code above will simply add an empty cell and go to the next line.
For Case 2 we can use unique to get a list of all define types, create an empty data.frame with those types and have a loop (or sapply) to add the values in the right columns.

We can check the nature of metadata.txt:

 types <- gsub("^(.*)\t.*$","\\1",dat)
 sapply(unique(types), function(t) length(grep(t,types)))

     STANDARD_NAME                ORGANISM            ORGAN_SYSTEM                    PMID       PUBLICATION_TITLE
                257                     257                     257                     257                     257
            AUTHORS                   GEOID            EXACT_SOURCE    EXTERNAL_DETAILS_URL                    CHIP
                257                     257                     257                     257                     257
      CATEGORY_CODE             CONTRIBUTOR         CONTRIBUTOR_ORG       DESCRIPTION_BRIEF RAW_PUBLICATION_MEMBERS
                257                     514                     257                     257                     257

Your are lucky and you don't have to worry about shifting because of missing values.
There is a fixable issue with CONTRIBUTORS, but it is related to the regex, not the file format.


which is the right prosite pattern? I would like to ask if the right Prosite pattern for this multiple alignment:

Is it this:

A-T-[AT]-G-x-C-[AGC]-C-x(1,4)-A

or this:

A-T-[AT]-G-x-C-[AGC]-C-x(1,3)-[ATGC]-A

Actually, I would like to ask if x always implies a gap, and if it is proper to use it when any symbol but not gap can be in the sequence. Thank you!

If nucleotides are acceptable in prosite format because I thought it was protein motifs, the most accurate would be,
A-T-[AT]-G-x-C-[AGC]-C-[GCA](0,1,3)-x-A

The key point is that,

[GCA](0,1,3) of x(0, 1, 3) means zero/ nucleotides (3 indels) or one/ G or C or A nucleotide, and 2 indels or three/ G or C or A nucleotides.
x is any nucleotide i.e. A or G or T or C - it doesn't mean an indel. x(0, 1) is the best representation of a single indel (you call them gaps), here meaning any nucleotide may or may not be present at a given site.

Example from Prosite documentation:

<A-x-[ST](2)-x(0,1)-V
This pattern, which must be in the N-terminal
of the sequence ('<'), is translated as:
Ala-any-[Ser or Thr]-[Ser or Thr]-(any or none)-Val


What is indel calling and what is its purpose? I'm having a difficulty in grasping the general purpose and concept of indel calling.

What exactly is this process?

There are two types of INDELs: short indels and long indels. Some put the threshold at 50bp; others choose 1000bp. Short and long indels are called differently.

<50bp short indels are called from read alignment directly. Modern indel callers essentially build a multi-alignment of reads and call an indel with enough read supports. Short indels may break open reading frames.

For Illumina reads, long indels are usually called with split alignment or read pair alignment because mainstream aligners are unable to directly align through a long indel. Long indels are a type of structural variation. They may affect gene structures in a more dramatic way (e.g. delete a whole exon; screw up transcription by transposon insertion; create pseudogenes).

Ultimately, both short and long indels are types of genetic variants. Calling them helps to understand how genetics shapes phenotypes.


Bioconductor Conda release delay How long does it take for the latest Bioconductor to end up on Conda? I think I saw a blog post about the challenges with this process, but I cannot find it. For example, Bioconductor 3.13 was released a week ago, but is still not available. For comparison, many CRAN packages tend to show up within a day.
Additionally, is there a way to keep track of that progress?

Devon will probably add a more precise answer, until then I just link his Tweet from last week:
https://twitter.com/dpryan79/status/1394714988720308226


No output when trying to obtain protein sequence from PDB file I am trying to obtain the sequence from PDB file, but no output is given when I run the following code :

from Bio import SeqIO

with open("1a9l.pdb", "r") as handle:
    for record in SeqIO.parse(handle, "pdb-atom"):
        print(record.seq)

That PDB entry is an RNA molecule, so there is no protein sequence to extract.

If it extracted the RNA sequence by default here you can see there could be an ambiguity as to whether it was an RNA sequence or a protein sequence, as the alphabets share characters.


How can I build a protein network pathway from a gene expression quantification file? Assume I have found the top 0.01% most frequent genes from a gene expression file.

Let's say, these are 10 genes and I want to study the protein protein interactions, the protein network and pathway.

I thought to use string-db or interactome, but I am not sure what would be a plausible way to approach this problem and to build the protein network etc. Are there other more suitable databases?

How can I build a mathematical graph or network for these data?

There are multiple ways to do this, and multiple protein interaction databases besides the ones you mentioned, such as BioGRID or IntAct. Interaction databases are different in how interactions are defined, sometimes it can be experimental evidence of interaction, sometimes coexpression, orthology-based predictions, etc.

There is no single solution to your problem. For String-DB you can use their R package STRINGdb.

# string-db
library(STRINGdb)

# retrieve full graph for H. sapiens
string_db <- STRINGdb$new(version="10", species=9606,
                          score_threshold=400, input_directory="" )

# define genes of interest
genes.of.interest <- ...

# get their neighbors
neighbors <- string_db$get_neighbors(genes.of.interest)

# get the subgraph of interest (your genes + their neighbors)
my.subgraph <- string_db$get_subnetwork(c(genes.of.interest, neighbors))

# look how many genes and interactions you have
my.subgraph

If you prefer GUI, you might consider using Cytoscape. They have multiple addons for interaction databases, e.g. for STRINGdb. There you can just provide your gene list and get a network in a few clicks.


Installation of PRANK MSA in WLS Ubuntu 20.04 LTS I want to install PRANK on the Windows 10 Linux subsystem (Ubuntu 20.04 LTS), I have followed the installation instructions to no avail.
wget http://wasabiapp.org/download/prank/prank.source.140603.tgz
tar xvzf prank.source.140603.tgz
cd prank-msa/src
make
./prank

[ cp prank ~/bin ]
[ sudo cp prank /usr/bin ]

I followed those instructions but got stuck with cd prank-msa/src I don't know why they want me to enter a directory that doesn't exist. I also followed this video on how to install linux files but it wasn't helpful. Any suggestions would be appreciated.

Their instructions are outdated and don't match what exists on their server. The problem you're running into is not at the cd step but even before that - at your first step even. Your wget fails because that file does not exist. Go to http://wasabiapp.org/download/ and see for yourself - there is no /prank/ directory in there.
I tried using curl, and it wrote the 404 page as a .tar.gz file, so if you did that, the tar step would error out.


Prevent Snakemake from changing wildcards: sample for every rule EDIT: I posted this as a bug, perhaps wait for the result of that? github

All,

I am using Snakemake to build a pipeline. To define input and output names, I constructed them based on a filename base (from the samples table), a read indicator (like R1, R2) and the fastq file extension, like so:

samples = pd.read_csv('samples.tsv', sep='\t').set_index('samples', drop=False)
reads = ['R1', 'R2']
fastq_file_extension = _001.fastq.gz

rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"

However, during pipeline execution already I can see strange things happening, for this rule it reports:

[Fri Feb 21 12:39:57 2020]
rule fastqc:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz
output: qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.html, qc/fastqc/0053_P2017BB3S20R_S2_R1_fastqc.zip
jobid: 10
wildcards: sample=0053_P2017BB3S20R_S2_R, read=1

As you can see, Snakemake pulled the R part of the read to the sample and set the read to 1 (should be R1). But this still works because it produces output files as Snakemake expects.

But the next rule, trim_galore_pe:

rule trim_galore_pe:
    input:
        [os.path.join(config['Directories']['fastq_raw_dir'], '{sample}' + read + fastq_file_extension) for read in reads]
    output:
        fastq_read1 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.1_val_1.fq.gz'),
        fastq_read2 = os.path.join(config['Directories']['fastq_trimmed_dir'], '{sample}.2_val_2.fq.gz'),
        report_read1 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.1.fastq.gz_trimming_report.txt'),
        report_read2 = os.path.join(config['Directories']['trim_galore_qc_dir'], '{sample}.2.fastq.gz_trimming_report.txt')
    params:
        extra="--retain_unpaired"
    wrapper:
        "0.49.0/bio/trim_galore/pe"

This goes wrong, during pipeline execution Snakemake reports:

[Fri Feb 21 12:40:37 2020]
Error in rule trim_galore_pe:
jobid: 7
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
conda-env: /home/nlv24077/projects/snaqs/snakemake_envs/1b94430b
cluster_jobid: 11904662.osios
Error executing rule trim_galore_pe on cluster (jobid: 7, external: 11904662.osios, jobscript: /home/nlv24077/temp/RNASeqPairedEnd/.snakemake/tmp.ooypuklc/snakejob.trim_galore_pe.7.sh).
For error details see the cluster log and the log files of the involved rule(s).

This is understandeble because the file fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz does not exist, it should be fastq_trimmed/0053_P2017BB3S20R_S2.1_val_1.fq.gz, somehow Snakemake added an extra underscore, which you can see in the execution output, under "wildcards" (0053_P2017BB3S20R_S2_).

[Fri Feb 21 12:39:57 2020]
rule trim_galore_pe:
input: fastq_raw/0053_P2017BB3S20R_S2_R1_001.fastq.gz, fastq_raw/0053_P2017BB3S20R_S2_R2_001.fastq.gz
output: fastq_trimmed/0053_P2017BB3S20R_S2_.1_val_1.fq.gz, fastq_trimmed/0053_P2017BB3S20R_S2_.2_val_2.fq.gz, qc/trim_galore/0053_P2017BB3S20R_S2_.1.fastq.gz_trimming_report.txt, qc/trim_galore/0053_P2017BB3S20R_S2_.2.fastq.gz_trimming_report..txt
jobid: 7
wildcards: sample=0053_P2017BB3S20R_S2_

So how can I prevent Snakemake from changing the samples and read wildcard in seemingly arbitrary ways?

Placing two wildcards directly behind each other will make it impossible for snakemake to determine where one wildcard ends and the other begins. So it is simply taking the longest possible match for the first wildcard (sample) and then one character for the second (read). This will also cause the _ to be included in the wildcard value sample.

I'm assuming that this same pattern is used for the inputs of another rule to refer to the trim_galore_pe output. Hence why the _ is included in the wildcards there.

You might be able to fix this, simply by adding a _ in-between {sample} and {readgroup} in your first rule (and any other rules using this same pattern):

rule fastqc:
    input:
        os.path.join(config['Directories']['fastq_raw_dir'], '{sample}_{read}' + fastq_file_extension)
    output:
        html=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.html'),
        zip=os.path.join(config['Directories']['fastqc_dir'], '{sample}{read}_fastqc.zip')
    wrapper:
        "0.49.0/bio/fastqc"


How does MinKNOW classify 1D reads as "pass" or "fail"? I have found a couple of sources1,2 that indicate that a read in a 1D² run is classified by MinKNOW as "pass" and put into the fastq_pass folder if both of the following conditions are met:

Both strands were read
The average Phred score of all bases in the read was >9

(Thus a read in a 1D² run goes into the fastq_fail folder if either its Phred is too low or if the second strand didn't follow the first through the pore.)

However, I have been unable to track down any source indicating how the pass/fail classification works for 1D runs (where only one strand is supposed to be read per read). One might reasonably guess that it's based solely upon whether the mean Phred score is above 9, but inspection of the data from a 1D run reveals that this is not the case. While all "pass" reads seem to have a mean Phred score of at least 9 (the lowest I see is 9.45), some fail reads have scores above 9 too (with the highest one being 12.9):

Thus it looks to me like this is not the condition used to classify 1D reads. But then what is?

The classification is indeed based purely on comparing a mean to a threshold... but it's not the mean of the Phred scores that represent base quality on a logarithmic scale, but rather the mean of the underlying quality values prior to being transformed into Phreds.

If I compute the "mean" quality for a read in this way - by converting each Phred in the read back to a probability, taking the mean of those probabilities, and then converting that mean probability to a Phred - I find that whether a read is a pass or fail is determined by comparing that score to a threshold of 7:

A mean quality score computed in this way is also included in the mean_qscore_template column of the sequencing summary file. If you open that file in a spreadsheet program (it's just a TSV) and sort by mean_qscore_template, you will find that the passes_filter column is always FALSE for mean_qscore_template values below 7 and always TRUE for values above 7.

(Hat tip to Wouter De Coster for successfully guessing that this was the answer in the comments and pointing me on the right path.)


Download proteomes from NCBI based only on binomial names I'm trying to download around 2000 proteomes from NCBI, but all I have is the genus and species e.g: lophodermium_seditiosum trichoderma_harzianum

Does anyone know a good way of doing this in script form?
My code at the moment is:

from Bio import Entrez
from Bio import SeqIO
import os

file_name = "proteome1.faa"
net_handle = Entrez.efetch(db="protein", id="trichoderma_harzianum", rettype="gp", retmode="text")
out_handle = open(file_name, "w")
out_handle.write(net_handle.read())
out_handle.close()
net_handle.close()
print("Saved")

for seq_record in SeqIO.parse("proteome1.faa", "genbank"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

But I think it wants a number as the id, I don't know if it's possible to convert species name to a numerical ID within the script, or if that's the best way of going about this.

The general procedure is explained here:

First step is to search the Taxonomy database with the organism name. Accepted common names usually work at all taxonomic levels. Use the scientific name or formal name if no results are obtained with the common name.

In code:

from Bio import Entrez

Entrez.email = 'user@example.org' # Put your email here

def get_tax_id(species):
    species = species.replace('_', '+').strip()
    search = Entrez.esearch(term=species, db='taxonomy', retmode='xml')
    record = Entrez.read(search)
    return record['IdList'][0]

if __name__ == '__main__':
    organisms = ['lophodermium_seditiosum', 'trichoderma_harzianum']
    taxids = [get_tax_id(organism) for organism in organisms]
    print(taxids)

This gives us the taxids for our species, respectively 128023and 5544. Now we can use thes to get the proteomes. For Lophodermium seditiosum we find two results, but for Trichoderma harzianum there are 65900 results, so you might want to tweak the query to filter the results. For example, if you only want the sequences appearing in PDB you add AND pdb[filter] so that the query becomes txid5544[Organism:exp] AND pdb[filter].

def get_proteomes(taxid):
    query = 'txid{}[Organism:exp]'.format(taxid)
    handle = Entrez.esearch(db='protein', term=query)
    result = Entrez.read(handle)
    if int(result['Count']) <= 2000:
        ids = ','.join(result['IdList'])
        proteomes = Entrez.efetch(db='protein', id=ids, rettype='fasta').read()
        return proteomes
    else:
        raise RuntimeError('Too many results!')

if __name__ == '__main__':
    print(get_proteomes('128023'))

Returns:

>ADI44294.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS

>ADI44293.1 actin, partial [Lophodermium seditiosum]
CKAGFAGDDAPRAVFPSIVGRPRHHGIMIGMGQKDS


Store and search on variant and annotation data I need to store a large number of variants with their annotations. I need to import dbSNP, as well as Clinvar, and ExAC etc. Also I have sample variant data from patients. I need to be able to search on genomic regions as well as annotations. And it needs to be scalable.
What would be a preferred data storage technology for such a case? What is the industry using? It needs to be accessible by APIs/REST.
We are using ElasticSearch, but that is proving to be expensive in terms of servers needed.

Something similar to this has been discussed in
The state, limitations and comparisons of large variant stores, but that post is a little old, I am sure there has been progress since then.

The question you linked is about storing genotypes. Your question is about keeping site information, which is different and much easier.

As to solutions, you can parse the information you want and store them in a standard SQL database. You can use UCSC binning for region queries. Modern SQL databases all provide R-tree indices, which can be adapted for region query and is probably better these days.


Separation of mixed plasmid DNA sequences post whole-plasmid sequencing Imagine a DNA sample containing a mixture of different intact plasmids. These samples are sequenced using either MiSeq or HiSeq sequencing. Would it possible to assemble these plasmids post-sequencing as they would have been when sequenced individually?

You would expect to have high coverage, given the plasmids are short, so de novo assembly would be likely very easy.  Given that each plasmid is present in different multiples, you would expect different coverage on each plasmid, so it might be best to approach it as a metagenome-type or transcriptome assembly, rather than a classic genome assembly.

Alternatively, Google search just told me there is a tool called plasmidSPADES designed for just your purpose.  SPADES itself is a reliable genome assembler, so I would imagine plasmidSPADES to be a good choice.  It does assume that there is also some whole bacterial genomes in the sequencing mix though.

As @terdon commented, it does depend on whether you have many identical copies of each plasmid, or many similar but not identical plasmids in the mix.


GATK documentation for required depth to reliably call heterozygous mutation in diploid organism? I'm looking for official GATK documentation (or a recent manuscript) that defines a general recommendation/requirement for sequencing depth to reliably call a heterozygous point mutation in a diploid organism (WGS). In this case, I'm working with humans. Specifically, suppose we're looking at position X in the genome and we want to classify (True/False) whether there was enough coverage to detect a heterozygous point mutation.

What coverage would that be? I'm curious about both: (1) UnifiedGenotyper and (2) HaplotypeCaller

I know general practice is to have 30x coverage genomewide, but I don't see anything in GATK's documentation (i.e., Best Practices) about required coverage.

Really appreciate your help.

Note: I also submitted a question to the GATK team end of last week.

UPDATE: As an update, Sarah Walker (co-author on the poster) responded to my question on the GATK forum. She clarified with the following statement:

  We believe the sites around 30X (and above 150X) are being filtered
  due to low mapping quality (since it is whole genome so there are many
  areas that are hard to map), which explains the low sensitivities in
  these areas. As the error bars show, there are very few variants in
  these low mapping quality areas.

She then posted the following plot that falls much more inline with what we would all expect. This is a plot showing sensitivity for NA12878 compared to NIST truth.

shlee, from the GATK forums, pointed me to a poster from AGBT 2018 where the authors compared the HiSeqX and NovaSeq instruments. Figure 3 shows sensitivity for SNPs and short INDELs across multiple allelic fractions (AF) from spike-ins. I was surprised to find that SNPs where AF = 0.4 only had ~50% sensitivity at ~25x coverage, and only ~75% sensitivity at ~50x coverage.

The poster is targeted at somatic variation, but data from figure 3 should behave like a standard diploid germline variant because you often see a wide range of allelic fractions from germline variants. They used spike-ins for figure 3. It's not clear to me whether they used HaplotypeCaller in figure 3, but I don't see how they could have used Mutect2, like they did in figures 1 & 2.

Based on this information, it appears the standard 30x genome-wide coverage is generally insufficient to confidently call many heterozygous mutations, and that ~90-100x is ideal.

Am I missing something? Are we really missing that many heterozygous mutations in standard WGS/WES?

I'll leave this question unanswered for now to allow others to chime in on interpretation, etc. I'd appreciate additional thoughts.


Is there a modern alignment tool tailored for transmembrane regions? I am looking for a project or tool that allows programmatic pairwise alignments of proteins but that takes care with transmembrane regions of proteins. TM regions are traditionally too information poor for many tools to align accurately (Wong et al, 2010).

TM-aligner doesn't appear to have a program to download, and is only available via webserver. On the otherhand,
PRALINE has some optimisation for TM regions, however was released in 2008 and since then there has been a boom in TM knowledge from structures so does not benefit from that.Currently I am using PSI-BLAST. However,
PSI-Blast uses SEG to mask low complexity regions. This will most likely mean that TMHs are masked making the alignment of those regions fairly unreliable.

Primarily I am interesting in sequence alignments. But if anyone knows of a TM structural alignment, that could be good to know.
I think generally that counting idea works for an ideal helix, but a vast majority of TMHs break convention. I wouldn't expect to see a pattern based aligner that was that general. I'd be interested if anyone knows of a tool that juggles that though!
If I understand correctly, Matteo was saying that TMHs follow a simple patter (HHP/H,HHP/H...~7 times). This is not the case. There are too many examples of TMHs that do not match this pattern to assume this pattern would accurately align TMHs. As @MichaelG. suggests, the number of TMP structures has not caught up to their soluble counterparts in the database. I've updated the question to explain why PSI-Blast might not be the best option.

Answer from @matteo-ferla converted from comment:
I assume you mean sequence (MSA), not structure alignment. Despite being a structural bioinformatician, I would slap a big caveat against structure based aligners. TCoffee:Espresso is probably the best structure-based aligner (browser and download and API), but the results are not as solid as some non-structural approaches such as training Psi-blast and aligning with that matrix. I am not familiar with MSA of TM regions, however, so this group may be an exception.
Given that helices are a repeat of 2 hydrophobic residues and 1.5 either polar (core-wards) or also hydrophobic (depending on how the helices pack), a pattern based aligner ought to do well from an evolutionary perspective. If you really want from a structure perspective and don't want to know about the cases where a loss of a turn at the top is compensated by the bottom, say, then TCoffee is definitely worth a try
Answer from @m converted from comment:
PSI-blast is a useful tool for deep alignment patterns, albeit it has its critics. I agree the vast majority of in vivo structures lack TM proteins, for different reasons


Seeking explanation of the hg38 files downloaded from bowtie 2 website I downloaded the H. sapiens, NCBI GRCh38 files from Bowtie's website. After unzipping, there are 6 files, 4 that end in set 1.ebwt, set2.ebwt, set3.ebwt, and set4.ebwt and two that end in set.rev.1.ebwt and set.rev.2.ebwt.

I am unable to find any explanation for what these files are either online.

which one do I use for bowtie2 mapping?

Thanks

Bowtie is providing you with the index files for hg38.  This are the result of the bowtie2-build indexer.  However you have downloaded the indexes for Bowtie1 not bowtie2

  Bowtie 2’s .bt2 index format is different from Bowtie 1’s .ebwt
  format, and they are not compatible with each other.

Then according to the manual you pass the basename of the index for the reference genome. So this should be the name prior to the extension, .1.bt2 for Bowtie2.

bowtie2 [options]* -x <bt2-idx> {-1 <m1> -2 <m2> | -U <r> | --interleaved <i>} -S [<sam>]
Main arguments
-x <bt2-idx>
The basename of the index for the reference genome. The basename is the name of any of the index files up to but not including the final .1.bt2 / .rev.1.bt2 / etc. bowtie2 looks for the specified index first in the current directory, then in the directory specified in the BOWTIE2_INDEXES environment variable.


Find pattern that is present twice and allow <=2 mismatches on each I have a fastq file of 400,000 reads (so speed is important). In the sequences there are barcodes integrated that should be present twice. Given a barcode, I want to find the sequences that have the barcode present twice with <= 2 mismatches. So, with a barcode 'ATTCGACCGATAGG', I would like to retrieve all of the following sequences-

TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA
TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA

Note that the first barcode in the fourth sequence is short of one character. I have tried with biopython and regex but it's just too slow given I have 5k barcodes. I am wondering if there is a fast solution available in python or in something like grep, awk or anything else. Thanks.

This isn't perfect but you can build on it to deal with any issues (ie overalpping hits). Also if you need it faster just chunk it and run in parallel threads.
I just copy pasted your examples seqs into files to read (called them seqs and code).
#need blast installed
import subprocess
import pandas as pd

path = '/Users/ur/dir/'
#make a blast db
with open (path + 'seqs.txt', 'r') as fin:
    with open (path + 'seqs.fa', 'w') as fout:
        for i, line in enumerate(fin):
            fout.write('>'+str(i)+'\n')
            fout.write(line)

subprocess.run([
    'makeblastdb',
    '-in', path + 'seqs.fa',
    '-dbtype', 'nucl'
])

#use blast
subprocess.run([
    'blastn',
    '-task', 'blastn-short',
    '-query', path + 'code.txt',
    '-db', path + 'seqs.fa',
    '-out', path + 'out.txt',
    '-outfmt', '6'
])

#grab seqs with 2 hits (ie ~2x len(query)
df = pd.read_csv(
    path + 'out.txt',
    sep='\t',
    header=None,
    names=['query', 'db_name', 'percent','len'],
    usecols = [0,1,2,3]
)
df = df.groupby('db_name').agg('sum')
print(*list(df.query('len > 27').index))

I this case they are all hits. If you look at the blast output you'll see one hit was fragmented (blast isn't perfect with small seqs). But total len works to resolve this.


Alignment based vs reference-free (transcriptome analysis)? I want to focus on transcriptome analysis. We know it's possible to analyze RNA-Seq experiment based on alignment or k-mers.

Possible alignment workflow:

Align sequence reads with TopHat2
Quantify the gene expression with Cufflinks

Possible reference-free workflow:

Quantify sequence reads with Kallisto reference-free index

Both strategy generate gene expression table.

Q: What are pros and cons for each of the approach? Can you give guideline?

First of all, I would emphasize that "alignment-free" quantification tools like Salmon and Kallisto are not reference-free. The basic difference between them and more traditional aligners is that they do not report a specific position (either in a genome or transcriptome) to which a read maps. However, their overall purpose is still to quantify the expression levels (or differences) of a known set of transcripts; hence, they require a reference (which could be arbitrarily defined).

The most important criterion for deciding which approach to use (and this is true of almost everything in genomics) is exactly what question you would like to answer. If you are primarily interested in quantifying and comparing expression of mature mRNA from known transcripts, then a transcriptome-based alignment may be fastest and best. However, you may miss potentially interesting features outside of those known transcripts, such as new isoforms, non-coding RNAs, or information about pre-mRNA levels, which can often be gleaned from intronic reads (see the EISA method).

This paper also has some good considerations about which tools may work best depending on the question you want to answer.

Finally, another fast and flexible aligner (which can be used with or without a reference transcriptome) is STAR.


about the normalization of RNAseq data for calculating distance for unsupervised learning I have been working on clustering using RNAseq data. To compute distance, what kinds of normalization is optimal?

Can we use normalization using relative log normalization? I understand this should be used for differential gene expression analysis, not for unsupervised learning in which we do not use information of disease condition.

If you can show good papers, please share.

Typically some sort of variance stabilizing transform is used before clustering. Popular options are regularized log transformation or a vst transform, which are available in DESeq2. Note that these are NOT used for performing differential expression, just things like clustering and PCA. For differential expression one would use TMM, RLE or similar.


How to get fasta alignment file from SAM/BAM file? I am not talking about consensus sequence, I know how to get consensus sequence using mpileup in samtools/bcftools. As I understand , SAM/BAM files are basically sequence alignment format so it's  natural to expect a straightforward way of conversion between them. Sadly I have yet to find one. Vcf-kit allows for conversion between vcf file and fasta alignment, but the fasta alignment contains only the SNPs. I need the whole alignment.  Does samtools/bcftools provide option for such conversion? If not, is there any other tool?

I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", you could mean pairwise alignment, but it is still impractical to output every pairwise alignment in a separate fasta file.

There are tools to convert SAM to BLAST-like format if that is what you want. The following shows one of them. I believe there are other tools easier to use, but I forget what are they.

Anyway, to install:

wget https://raw.githubusercontent.com/lh3/minimap2/master/misc/paftools.js
curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8

Then you can convert your ALN.sam with:

./k8 paftools.js sam2paf -L ALN.sam | ./k8 paftools.js view -l0 - | less -S

Or generate a SAM stream from ALN.bam and pipe it with:

samtools view -h ALN.bam | ./k8 paftools.js sam2paf -L - | ./k8 paftools.js view -l0 -

Here is an example output:

>MT_orang   16499   0   16025   +   MT_human    16569   576 16569   13773   16095   60  tp:A:P  mm:i:2150   gn:i:172    go:i:101
Ref+:        577 GTTTATGTAGCTTACCTCCT---CAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCA
                 ||||||||||||||  | ||   ||||||||| |||||||||||| | ||||||| |||| | |||||||||||||||||||||||||||||||||||||||||||||||| || |||||||||||||||||||||  |||||||| || ||||| || ||||   ||  || |||| |||||||||
Qry+:          1 GTTTATGTAGCTTA--TTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACA-CGCCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACATGCAAGCATCCCCGCCCCAGTGAG-TCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCATCA
//

You can change -l0 to -l80 to print alignments in multiple lines.


Understanding the initialization of the Needleman-Wunsch algorithm I would like to ask why when we initialize 0 as the starting point, there are 3 gray boxes that are not used as shown (highlighted in the red box) in the pic?

 Are these boxes to show that they are gaps or just to align them properly?

The first column is added in order to be able to align the sequences. They might be gaps if the alignment ends up beginning in different positions than the first. Without the first empty row and column it wouldn't be possible to have gaps at the beginning of the alignment

The dynamic programming is useful because you need first to calculate the score for each position, and then once you have the score, you can calculate the optimum path. So you solve your problem by first calculating what would happen if (there was a gap, there was a mismatch...) and then you find the optimum alignment.


Splitting fasta file into smaller files based on header pattern I have to split this fasta files into smaller files and write them into individual files my files

>lcl|CP000522.1_prot_ABO13860.1_1 [locus_tag=A1S_3471] [protein=hypothetical protein] [protein_id=ABO13860.1] [location=1..957] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13850.1_2 [locus_tag=A1S_3461] [protein=DNA replication protein] [protein_id=ABO13850.1] [location=950..1504] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13851.1_3 [locus_tag=A1S_3462] [protein=hypothetical protein] [protein_id=ABO13851.1] [location=complement(2523..3437)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13852.1_4 [locus_tag=A1S_3463] [protein=YPPCP.09C-like protein] [protein_id=ABO13852.1] [location=3538..4788] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13853.1_5 [locus_tag=A1S_3464] [protein=Cro-like protein] [protein_id=ABO13853.1] [location=5039..5629] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13854.1_6 [locus_tag=A1S_3465] [protein=hypothetical protein] [protein_id=ABO13854.1] [location=complement(6340..6906)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13855.1_7 [locus_tag=A1S_3466] [protein=Resolvase] [protein_id=ABO13855.1] [location=complement(7074..7685)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13856.1_8 [locus_tag=A1S_3467] [protein=hypothetical protein] [protein_id=ABO13856.1] [location=complement(8602..9732)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13857.1_9 [locus_tag=A1S_3468] [protein=putative lipoprotein] [protein_id=ABO13857.1] [location=complement(10072..10374)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13858.1_10 [locus_tag=A1S_3469] [protein=Diaminopimelate decarboxylase] [protein_id=ABO13858.1] [location=complement(10367..10723)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13859.1_11 [locus_tag=A1S_3470] [protein=regulatory protein LysR] [protein_id=ABO13859.1] [location=complement(12076..12444)] [gbkey=CDS]

The other pattern is

>lcl|CP000523.1_prot_ABO13861.1_1 [locus_tag=A1S_3472] [protein=DNA replication protein] [protein_id=ABO13861.1] [location=1..951] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13862.1_2 [locus_tag=A1S_3473] [protein=hypothetical protein] [protein_id=ABO13862.1] [location=3048..4262] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13863.1_3 [locus_tag=A1S_3474] [protein=hypothetical protein] [protein_id=ABO13863.1] [location=4357..5133] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13864.1_4 [locus_tag=A1S_3475] [protein=hypothetical protein] [protein_id=ABO13864.1] [location=6197..8608] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13865.1_5 [locus_tag=A1S_3476] [protein=secretory lipase] [protein_id=ABO13865.1] [location=8705..9403] [gbkey=CDS]

So now my idea is how do i parse and write them into individual files such that CP000522 output written to one file and CP000523 written to another file so forth and so on.

So far what i understand is i have to match the pattern after >lcl
so there are other patterns like "LN997847" in the file

Not sure how to proceed tried it in R but failed

it can be done with awk and sed which i tried but i can;t define something that parse all header like takes into account CP as well as LN .

Any suggestion or help would be highly appreciated .
here is my file


Principle of TMT Tags in Multiplex Proteomics I am new to proteomics research and analyzing mass-spec data. I am trying to learn more about the theory/principle of TMT tagging (a type of isobaric labeling) coupled with LC-MS/MS spectrometry. I understand that each tag has a

1) protein binding group (amine reactive ONLY?)
2) reporter group (so we know the tag came from which sample?)
3) normalisation group (to ensure that each tag has same mass?)
4) cleavage group (so tag can be cleaved from the peptide?)

Did I get the purpose of all the four parts in the tag?

Also, what does it mean by 4-plex, 10-plex, n-plex tag? And why do the tags have different isotopes?

If I could have some clarity on the workflow of TMT tagging and how it leads to analysis, it would be really helpful.

TMT stands for tandem-mass-tag. The idea is to tag different samples with one of the tags then pool all samples together and run them through the mass spectrometer together. Then the reporter groups are split off and quantified (ideally) for each peptide individually. There are different methods with different numbers of isobaric labels. Isobaric means all the labels have the same mass therefore they behave exactly the same in the mass spectrometer at first. E.g. TMT11 means you have 11 different tags with identical masses, so you can pool and process 11 samples at a time.
The protein binding group attaches the label to the peptide. Then the tagged peptides from all samples are run through the chromatography and the spectrometer. At first, the MS instruments simply collects the tagged peptides with different m/z ratios at each point in time for a while. When enough peptides are collected reporter groups are cut off and the peptides then quantified. This cutting-off happens by increasing the internal energy of the collected peptides. Because the cleavage group has the weakest chemical bond it brakes first and releases the reporter group. Again, this step separates the reporter group from the rest. Finally, the cut-off reporter groups are quantified.
The different reporter groups don't have the same mass, otherwise you would not be able to differentiate between them in the MS instrument. But you need same masses before quantification to ensure that the same peptides are all collected together. Therefore, you need to balance out the different masses of the reporter groups. That is the job of the normalization group. So, the normalization group ensures that overall the tags have the same mass and the peptides with different tags are not treated differently in the MS instrument before quantification. Then for quantification the reporter groups (with different masses between each other) are cut off and the MS instrument can therefore detect and quantify them individually. There is a time dependency of the whole process which happens in multiple phases. The whole process makes it possible to compare the amounts of peptides coming from the different samples.
The chemistry in the chromatography column is quite complicated which makes it difficult to compare samples from different experiments. The presence of one peptide might enhance or suppress the signal of another peptide or impact the chromatographic retention time. Therefore, you should only compare within the same sample. More tags TMT11, TMT16 etc. allow you to run more samples in one run.
You can then use one (or more) channels (mass-tags), to tag a standard sample and then compare everything to this standard that is assumed to be identical over multiple experiments. This way you can compare even more samples than available tags. e.g. more than 11.


Vegetable databases usage Good morning,
I'm aligning with blast some species, against all nt database to analyse the vegetables. However, in some cases we have species of animals mixed with vegetables and want only to detect the vegetable species. Normally, we would only expect to detect the vegetables only. But this is normal since we are using the all nt. So now, could anyone point me out, the way to detect only vegetables? Thanks for the help.

Using the blast command line tool, you can limit your search to a specific taxon using the –taxids option.

So, in your case -taxids 3193 should limit your search to plants.

If you are using the ncbi blast webinterface, you can restrict your search by entering "plants" into the Organism field in the "Choose Search Set" section


What is 'k' in sequencing? When a DNA sequence is sequenced, I've only ever dealt with A,T,C,G and N which indicates un-identifiable bases. However, I came across a 'k' recently and I had asked another researcher who gave me an answer for what 'k' represents but I don't quite recall. It can't just be an anomaly in that one sequence file. Any ideas?

See IUPAC codes:

So, as you can see above, K means "Either G or T".


genes with highest RPKM/FPKM for a human RNA-seq experiment? Which genes/transcripts show up as having the highest RPKM/FPKM on a human RNA-seq experiment? Is it usually the same genes/transcripts?

EDIT: mRNAs (e.g. poly-A RNA-seq preps) and not including small RNAs.

As asked as it is, the answer is probably no.
Indeed, the highest expressed RPKM/FPKMs will be different from one condition (or tissue) to another one for example.
Then, you may also have technical artefacts due to the wet-lab part or to the normalisation. For example, mitochondrial genes are often reported in the top expressed genes.
Now, if you want to compare your results to an outside source, I recommend checking the EBI gene expression Atlas* and more specifically the Baseline expression experiments.
Pro-tips:

You can filter the genes based on their nature (eg, mRNA) and/or their expression.
You can also visualise the overall expression of a specific gene across several experiments.

*Disclaimer: I am related to this project myself.