Profile of conservation between bacterial sequences and human protein My goal is to make a conservation plot between bacterial sequences and a human protein.

So far, I have a FASTA file of the protein, and a FASTA file with the sequences of the proteins from the BLAST results.

In my initial attempts, I have tried to do this using msa software:

I have downloaded clustalx and run a profile alignment using the single protein sequence as Profile 1 and the sequences from the BLAST results as Profile 2, and selected the option "Align Sequences to Profile 1".

I'm having difficulty interpreting and dealing with the results. Since the majority of the sequences aren't aligned, the results are a total mess, with a huge number of dashes added.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

ie:

I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.

Is there a better way to achieve this or to salvage the results from the profile alignment?

Thanks.

I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the same genome so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.

You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

  My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.

If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py

Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.

For example, given this MSA:

    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
           ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
           ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
           EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
           DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH

You'd get this plot:

  I would like to have a graph with my protein on the x-axis and a plot
  like the regions of high conservation for the multiple alignments,
  except with the spikes corresponding to a high number of BLAST hits.
  This would allow me to identify regions of my protein that have higher
  sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.

What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.

As you say in the comments,

  My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

Blast sequence of interest.
Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
Align all the sequences.
Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.


How to compute the Shannon entropy for a strand of DNA? I'm confused by the computation of sequence logo. Wikipedia gives a process about this without a concrete example.

Let's just consider DNA, so there are 4 bases (nucleic acids).

The following data comes from the book "Machine Learning - A Probabilistic Perspective (Figure 1)"

The Shannon entropy of position $i$ is:

${\displaystyle H_{i}=-\sum _{b=a}^{t}f_{b,i}\times \log _{2}f_{b,i}}$

Where ${\displaystyle f_{b,i}}$ is the relative frequency of base

This post is computing position 3, where it seems that the relative frequency of a is 100%, and the relative frequency of 3 other bases are 0%.

Now, how do you calculate the Shannon entropy of position 3? It seems to be 0, which is apparently incorrect, where am I wrong?

Why do you think the entropy of 0 is incorrect? It intuitively makes sense, as there is no uncertainty about the base at position 3, and thus there is no entropy.

However, what is plotted in a sequence logo isn't the entropy, but rather a measure of the "decrease in uncertainty" as the sequence is aligned. This is calculated by taking the entropy at this position if we randomly aligned sequences ($H_g(L)$), and subtracting from it the entropy of the alignment ($H_g(s)$):
$$
R_{sequence}=H_g(L) - H_s(L)
$$

The entropy at position 3 based on a random alignment is calculated by assuming there are 4 equally likely events (one for each base), and thus:
$$
\begin{align*}
H_g(L) & = -((1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2)) \\
H_g(L) & = 2
\end{align*}
$$
Notably, $H_g(L)$ will always be 2 when dealing with nucleotide sequences.

So in your example, we have:
$$
\begin{align*}
R_{sequence}&=H_g(L) - H_s(L) \\
R_{sequence}&=2 - 0 \\
R_{sequence}&=2 \text{ bits of entropy}
\end{align*}
$$

Finally, to work out the height of each base at each position in the logo, we multiply the frequency of that base by the overall information gain $R_{sequence}$. In this case the frequency of the base A at position 3 is of course 1:
$$
\begin {align*}
H(b,l) &=f(b, l) \times R_{sequence} \\
H(A, 3) &= f(a, 3) \times 2 \\
&= 1 \times 2 \\
&= 2
\end{align*}
$$

Source: Sequence logos: a new way to display consensus sequences.
 (Thomas D. Schneider, R. Michael Stephens, 1990)

Note: I've excluded the correction factor $e(n)$ from the overall information-gain calculation for the sake of simplicity. Refer to the paper above for an explanation of this.


Random access on a FASTQ file I would like to select a random record from a large set of n unaligned sequencing reads in log(n) time complexity (big O notation) or less. A record is defined as the equivalent of four lines in FASTQ format. The records do not fit in RAM and would need to be stored on disk. Ideally, I would like to store the reads in a compressed format.

I would prefer a solution that does not require any extra files such as for example a reference genome.

The title of this question mentions a FASTQ only because FASTQ is a common format for storing unaligned reads on disk. I am happy with answers that require a single limited transformation of the data to another file format in time complexity order n.

Update

A clarification: I want the random record to be selected with probability 1/n.

Arbitrary record access in constant time
To get a random record in constant time, it is sufficient to get an arbitrary record in constant time.
I have two solutions here: One with tabix and one with grabix. I think the grabix solution is more elegant, but I am keeping the tabix solution below because tabix is a more mature tool than grabix.
Thanks to user172818 for suggesting grabix.
Update
This answer previously stated that tabix and grabix perform lookups in log(n) time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of n in complexity. However, both tools use an index that scales in size proportionally to n. So, the loading of the index is order n. However, if we consider the loading of the index as "...a single limited transformation of the data to another file format...", then I think this answer is still a valid one. If more than one record is to be retrieved, then the index needs to be stored in memory, perhaps with a framework such as pysam or htslib.
Using grabix

Compress with bgzip.
Index the file and perform lookups with grabix

In bash:
gzip -dc input.fastq.gz | bgzip -c > output.fastq.gz

grabix index output.fastq.gz

# retrieve 5-th record (1-based) in log(n) time
# requires some math to convert indices (4*4 + 1, 4*4 + 4) = (17, 20)
grabix grab output.fastq.gz 17 20

# Count the number of records for part two of this question
export N_LINES=$(gzip -dc output.fastq.gz | wc -l)

Using tabix
The tabix code is more complicated and relies on the iffy assumption that \t is an acceptable character for replacement of \n in a FASTQ record.  If you are happy with a file format that is close to but not exactly FASTQ, then you could do the following:

Paste each record into a single line.
Add a dummy chromosome and line number as the first and second column.
Compress with bgzip.
Index the file and perform lookups with tabix

Note that we need to remove leading spaces introduced by nl and we need to introduce a dummy chromosome column to keep tabix happy:
gzip -dc input.fastq.gz | paste - - - - | nl | sed 's/^ *//' | sed 's/^/dummy\t/' | bgzip -c > output.fastq.gz
tabix -s 1 -b 2 -e 2 output.fastq.gz

# now retrieve the 5th record (1-based) in log(n) time
tabix output.fastq.gz dummy:5-5

# This command will retrieve the 5th record and convert it record back into FASTQ format
tabix output.fastq.gz dummy:5-5 | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'

# Count the number of records for part two of this question
export N_RECORDS=$(gzip -dc output.fastq.gz | wc -l)

Random record in constant time
Now that we have a way of retrieving an arbitrary record in log(n) time, retrieving a random record is simply a matter of getting a good random number generator and sampling. Here is some example code to do this in python:
Using grabix
# random_read.py
import os
import random

n_records = int(os.environ["N_LINES"]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system("grabix grab output.fastq.gz {0} {1}".format(rand_record_start, rand_record_end))

Using tabix
# random_read.py
import os
import random

n_records = int(os.environ["N_RECORDS"])
rand_record_index = random.randrange(0, n_records) + 1
# super ugly, but works...
os.system(
    "tabix output.fastq.gz dummy:{0}-{0} | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'".format(
        rand_record_index)
)

And this works for me:
python3.5 random_read.py

Disclaimer
Please note that os.system calls a system shell and is vulnerable to shell injection vulnerabilities. If you are writing production code, then you probably want to take extra precautions.
Thanks to Chris_Rands for raising this issue.


Why does `pratchet()` in `phangorn` not find the most parsimonious trees? I called pratchet() on a dataset with binary values for the characters and it returned a parsimony value of 35. I then calculated the parsimony value for the dataset with a tree that I myself constructed and it was 7. How is this possible at all? I thought that pratchet() was supposed to return the most parsimonious tree on the basis of the data set, but that clearly did not happen in this case. Is there other software that is guaranteed to find the most parsimonious tree?

Perhaps the tree which you have calculated is not being caught in the search heuristic? Have you tried (as a trial) setting the tree which you calculated as the starting topology, just to check that it is indeed being calculated on by pratchet()?

Alternatively, you could try playing around with the search heuristic settings, like varying from ratchet to NNI, or you could also try using an external program, like PAUP*.


Seurat FindMarkers() output interpretation I am using FindMarkers() between 2 groups of cells, my results are listed but i'm having hard time in choosing the right markers. Do I choose according to both the p-values or just one of them? If one of them is good enough, which one should I prefer? I am interested in the marker-genes that are differentiating the groups, so what are the parameters i should look for?
             p_val avg_logFC pct.1 pct.2  p_val_adj
UBD   3.080608e-06 1.7755753     1 0.000 0.07683959
CFB   1.262305e-05 1.3989233     1 0.067 0.31485675
RGS13 1.548593e-05 0.9009480     1 0.200 0.38626551
PYDC1 1.896309e-05 0.9622537     1 0.133 0.47299636
BIRC3 4.099414e-05 1.0129472     1 1.000 1.00000000
ICAM1 4.379895e-05 0.8219610     1 0.200 1.00000000

as you can see, p-value seems significant, however the adjusted p-value is not. and when i performed the test i got this warning In wilcox.test.default(x = c(BC03LN_05 = 0.249819542916203,  ... : cannot compute exact p-value with ties
I am completely new to this field, and more importantly to mathematics. Please help me understand in an easy way.

The p-values are not very very significant, so the adj. p-value. You need to plot the gene counts and see why it is the case. It could be because they are captured/expressed only in very very few cells. VlnPlot or FeaturePlot functions should help. You haven't shown the TSNE/UMAP plots of the two clusters, so its hard to comment more.


How to download multiple proteomes at once? I'm looking for a way to download multiple proteomes at once from one clade, as fasta. Probably from the NCBI, because it looks more user-friendly than Ensembl, but Ensembl is okay too. In the best case scenario, I would like to constrain the download to a specific clade (eg. Actinobacteria, Chloroflexi etc.), maybe with a specific RefSeq Category (eg. reference/representative genome) or even depending on the genome assembly level. Command line would probably be the best option, but if there's a direct download option it would be nice to know it as well.
Does anything like this exist? Help would be much appreciated!

You could try download-refseq-genomes. It will fetch all genomes from the NCBI FTP server that are in a specific subtree of the phylogeny.

For example, downloading the amino acid sequences from all proteins in each genome assembly for all species belonging to Chloroflexi (NCBI taxon ID = 200795):

download-refseq-genomes.pl -t faa 200795

This will download all assemblies with status=Complete. Use option -a to download all assemblies with any status.


no result from heat map WGCNA I have created a heat map using WGCNA using R. But the heat map command is not showing the plot and there is no error. I have a large number of modules and traits. I am not sure if that is the reason. But I tried making it in HPC R and it still is not showing any result plot. Thank you!

The plot is probably being drawn to a device that you cannot view on a cluster. If you're using the plotNetworkHeatmap function, say, try something like this:
png("plot.png", ...) #specify height, width, etc
plotNetworkHeatmap(....) #call your function
dev.off()

Now, scp the plot.png file over to your local machine to view the plot.


What value of reference RMSD after docking simulation from any starting coordinate would be considered outrageous to have? I need to lightly verify a set of docking simulations from the reference RMSD of the outputs, and as we all know reference RMSD depends on the initial coordinates of the ligand set. What value of reference RMSD after docking from any starting coordinate would be considered outrageous to have?

Let's get some definitions out of the way:

RMSD, root mean square deviation, is a metric of distance of molecule A and molecule B. Think of it as average euclidean distance. It us casually called just 'deviation', but has nothing to do with standard deviation (the root of variance, 2nd moment). The units are Ångström.

Reference RMSD is a term that appears in some manuals and is a the RMSD between the calculated pose and the starting pose. When you set up your docking you as a user define the starting position and how much it can deviate from that position. If you have a crystal structure with a ligand and you are docking the same or similar compound, then yes, it makes total sense to worry about this metric. In which case, RMSD greater than 2 are worrisome in a 2 Å resolution structure. You especially see RMSD as a metric of accuracy when PDB structures are redocked, wherein the same ligand is docked and reference is the file from the PDB (or its PDB Redo derivative), so RMSD greater than 2 Å is bad as it would not match the density, while under 1 is really good —note that in tests about 50% and 75% of re-docked ligands are docked less than 2 with AutoDock (free) and ICM-Dock (pay) respectively. But if you are doing it for an experiment and not for a test, you would probably have set the docking to be highly constrained in the first place (the grid size in AutoDock or position randomisation factor in most others, eg. Rosetta Ligand dock). If you actually care about reference RMSD the thing to manually look at is the top 3-4 poses in the electron density of your native structure. If you don't really know how your ligand docks please ignore this value as it is just a product of your starting pose. It is not a measure of how "happy" the complex is.

ΔGibbs is a negative value in kcal/mol that is basically a metric of how "happy" the atoms are. It is actually an energy potential, but confer Wikipedia pages for the proper definition. The more negative it is the happier the atoms are. Force-field calculations solve this value and the difference between undocked and docked (ΔΔG) is a good metric of how well your ligand docks. Please use this.

Affinity ΔΔG is not a k_off or a k_on. That is a result of the ligand binding to nearby positions and how well it can get into the binding site (e.g. in P450 enzymes the ligand orientation is dictated by how it squeezes into the entrance way). And how insoluble the ligand is, which is a major problem*. For a better metric, there is dynamic undocking where a constrained MD simulation is done with the ligand getting pulled out of the pocket. But this is not trivial.

∗ As any comp-chemist can tell you it is mighty hard to calculate solubility. Some force-field calculations do not even take account of π-π or S-π interactions. Drug-like compounds are rigid aromatic and rigid and highly hydrophobic. Many programs use a hack called implicit waters, where the water is not simulated as individual molecules but as a flowing thing. This means that ΔΔG benefits of constrained waters are ignored and alternative protontation states of water too. This is fine for protein (mostly), but has major errors for drugs. But in terms of thing to look out for in a docking algorithm, explicit waters is not the first as conformer sampling and side chain perturbation are more important.


Entrez.efetch returns incomplete genbank records I am using the biopython Entrez.efetch command to retrieve all features (CDS, mRNA, ...) of some genomes.

In this case (NC_014649, Acanthamoeba polyphaga mimivirus), it works as expected:

from Bio import Entrez, SeqIO
handle = Entrez.efetch(id='NC_014649', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features

[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1181549), strand=1), type='source'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='gene'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='mRNA'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2579), strand=1), type='CDS'),
...

However in this case (NC_010830, Candidatus Amoebophilus asiaticus), it only returns one record for the whole genome.

handle = Entrez.efetch(id='NC_010830', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features
[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1884364), strand=1), type='source')]

However, if I check the corresponding file in the browser, it does have features: https://www.ncbi.nlm.nih.gov/nuccore/NC_010830.1

What could be the reason for this inconsistency ? Did I  miss something ?

This has been answered before. See: https://stackoverflow.com/a/55402322/6262370

In short, you need to either use rettype='gbwithparts' or rettype='gb', style='withparts' to download the entire genbank flat file.


How can I extract all known mutations of my BAM (or SNP/INDL files)? I am using a Genome Explorer tool to see all the mutations on my own DNA.

My particular interest is on listing the variants and get their names/ids.

Here are a few screenshots:

You can see that they tell you the genomic position, the allele frequency of the mutation, the name of the disease, the Category, etc etc...

However, the filters on this tool work very badly so I cannot filter all my mutations only...

So I was trying to find out an alternative tool to get this data in such a nice way this tool provides. I tried IGV but they just provide raw readings, I didnt find a way to get a list of variants the same way this Genome Explorer does.

The closest thing I've found was Geneious:

but it crashes on my computer and is too complex to use (I wasn't able to reach that screen above, from my BAM file). And it also eats a lot of RAM that I don't have.

Are there any alternatives? Even paying?

If you have a BAM file, you could easily run a variant caller analysis software such as freebayes, though you would need a good amount of computing power like a cluster or a dedicated node with good memory. Freebayes outputs a vcf file that has all the variants listed. Then you can load it into the UCSC or Ensembl human genome browsers and look at what genes they are.


How to generate a weblogo image that takes into account sequence frequencies? I have a large aligned amino acid sequences file (fasta file) from next generation sequencing data.
The unique sequences have been aligned but these all have certain read counts (frequencies associated with them).

I would like to generate a weblogo image that has the true probability of amino acids at each position.
Is there a way to input frequencies/read counts with the aligned files of sequences using weblogo?
For example, in the title of the fasta file or as a csv/tsv file?

So that when generating the weblogo image it takes into account the abundances of each sequence?

Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.

Generate a simple mapping file (counts) with the number of times each sequence was found:

seq1 3
seq2 10
seq3 2

Then, given an aligned protein fasta file (foo.aln) like this:

>seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

You can simply repeat the sequences as needed (see my answer here for the FastaToTbl and TblToFasta commands):

$ awk '
        (NR==FNR){
            count[$1]=$2;
            next
        }
        {
            for(k=1;k<=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts <(FastaToTbl foo.aln) | TblToFasta
>seq1.1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.2
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.3
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2.1
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.3
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.4
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.5
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.6
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.7
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.8
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.9
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.10
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3.1
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
>seq3.2
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

Note that the above assumes your shell supports the <(command) syntax for process substitution. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run FastaToTbl foo.aln > foo.tbl first, and then use foo.tbl as input.


Running python script from terminal: umap issue Some time ago I asked the question:

RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'

It was solved by following the answer in it, however, it is still not working when I run a python script from terminal, giving the same error:

python python/umap.py -path
"/Users/nikitavlasenko/python_scripts/output/" -data
"/Users/nikitavlasenko/python_scripts/data_files/WT_SCID_CCA.csv"

  Traceback (most recent call last):
    File "python/umap.py", line 26, in
      result = umap.UMAP().fit_transform(data.T)
  AttributeError: module 'umap' has no attribute 'UMAP'

From inside spyder the same script runs perfectly, but somehow when running it from the terminal it is not successful. I am new to python, so maybe I am not pointing correctly to umap library? Maybe it should be imported in a different way? Here is the whole script:

import umap
import numpy as np

def getopts(argv):
opts = {}  # Empty dictionary to store key-value pairs.
while argv:  # While there are arguments left to parse...
    if argv[0][0] == '-':  # Found a "-name value" pair.
        opts[argv[0]] = argv[1]  # Add key and value to the dictionary.
    argv = argv[1:]
return opts

if __name__ == '__main__':
  from sys import argv
  myargs = getopts(argv)
  print(myargs)
  path_to_save = myargs['-path']
  path_to_data = myargs['-data']
  data = np.loadtxt(open(path_to_data, "rb"), delimiter=",", skiprows=1)
  result = umap.UMAP().fit_transform(data.T)
  np.savetxt(path_to_save + "/umap.csv", result, delimiter=",")

Python version that is used by spyder and in the terminal seems to be the same: 3.6. I am using Mac OSX.

You answered your own question, yet it is not a coincidence that you have to change the name of your script.

When Python executes your custom script, it reads:

  import umap

Python will first search the local sources (the directory of the script) , and afterwards the installed sources.

In this case Python identifies the local source to be umap.py, and imports the (currently executed) script, while foregoing the installed sources.
Since you are calling functions, not defined in your current script, Python will not be able to fully execute your script.

For more information on importing see Module search path at docs.pyhton.org


efetch nucleotide -> protein ids I'm downloading translated viral genomes initially via Blast - which had shortcomings - and now by efetch.

What I want to do is obtain a complete list of all translated protein ids from a nucleotide Genbank entry.

efetch -db nucleotide -format gpc -id MT039887 | xtract -insd ....

The protein id and the translated gene sequences are within the efetch -db nucleotide -format gpc -id MT039887 output ... but how do I use xtract to parse them?

If you want the cds in fasta format, try:

esearch -db nucleotide -query MT039887 | efetch -format fasta_cds_na

And for the proteins sequence in fasta format:

elink -db nucleotide -target protein -id MT039887 | efetch -format fasta


Perfect Phylogeny vs Maximum parsimony I am searching various sources about phylogenetics. I saw some materials about perfect phylogeny and also phylogenies acquired from maximum parsimony constraint. They seem very similar to me. Are they the same?

No. A perfect phylogeny is one such that all characters evolve on the tree with no homoplasy, i.e. for a binary character, changes occur once from 0 -> 1 but never from 1 -> 0. A maximum parsimony phylogeny may produce a perfect phylogeny, but typically for real datasets some degree of homoplasy is required to explain character patterns.


Get frequency of index sequences of fastq file I recently did some sequencing on a MiniSeq and unexpectedly the Undetermined_S0_L001_R1_001.fastq.gz file contained quite a lot of data (about 40% of total). I would like to know which indices are present in this file and ideally get a frequency of all index sequences in this file. How can I do this?
This is the output of head:
@MN00148:128:000H2YCW2:1:11102:12116:1065 1:N:0:0
GNAGAGGATTCAGGGAGCCAGCGCGCATATATCAG
+
F#FFFFFFFFFFFFFFFFFFA=/A///////////
@MN00148:128:000H2YCW2:1:11102:11354:1065 1:N:0:0
GNACCCGTGCGGCTGGAGAGACCACTAGAAAGAAT
+
A#AFFFFFFFFFFFFFFFFF///////FF=////6
@MN00148:128:000H2YCW2:1:11102:4103:1066 1:N:0:0
ANGTTGGTGGTACATCAGGGACGCGGAGATCTCAG

I actually went back to the output folder on the MiniSeq and found this file: DemultiplexSummaryF1L1, which contains demultiplexing information of the whole run. Half-way through this file I found what I was looking for:
### Most Popular Index Sequences

### Columns: Sequence ReverseComplement HitCount

Index
TGACCA  TGGTCA  2455964
CGATGT  ACATCG  2320579
GGCCCC  GGGGCC  2023652
CGCGGG  CCCGCG  1983605
CCGTCC  GGACGG  1641274
CAGATC  GATCTG  1514252
CCGGCC  GGCCGG  1319383
CCGCGC  GCGCGG  1109940
ACAGTG  CACTGT  891344
CCCGGG  CCCGGG  706083
TGCCCA  TGGGCA  168294
CAGCTC  GAGCTG  164091
.... [data continues]


Transform and feed data into DESeq2 with DESeqDataSetFromMatrix There is a normalized expression matrix. I split it into two and want to do DE on the two cells' subsets. I am having trouble transforming it into the format that DESeq2 would accept. The function that I would think I need to use is the following:

dds <- DESeqDataSetFromMatrix(countData = cts,
                          colData = coldata,
                          design= ~ batch + condition)

It would be perfect if I could somehow feed two matrixes separately, but I can not find it. Any suggestions would be greatly appreciated.

  Update

I tried the following:

cells <- c(cell_names, cell_names_cyt)
groups <- c(rep("NE", times=length(cell_names)), rep("EX",
times=length(cell_names_cyt)))

colData <- data.frame(
  Cells <- cells,
  Groups <- groups
)

dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

But getting the error:

  Error in DESeqDataSet(se, design = design, ignoreRank) :
    some values in assay are not integers

Of course, they are floats since I already normalized the matrices.

  Update

I did the following since I am using Seurat:

df <- as.data.frame(as.matrix(seurat_object@raw.data))
df <- df[,colnames(df) %in% cells]

But then the following:

dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

Gave me an error:

  Error in DESeqDataSet(se, design = design, ignoreRank) : all variables in design formula must be columns in colData

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables
  in design formula must be columns in colData

I think that means just what it says...you don't have a "Groups" column in coldata.  Did you spotcheck your ColData data frame to make sure it looks the way it should?

> cells <- c("1","2","3")
> groups <- c("a", "b", "c")
> test <- data.frame(Cells <- cells, Groups <- groups)
> test
  Cells....cells Groups....groups
1              1                a
2              2                b
3              3                c
> colnames(test)
[1] "Cells....cells"   "Groups....groups"


Clustering using cellranger Is it possible to run cell ranger clustering using only a gene x cell matrix? Cellranger reanalyze performs clustering but requires a hd5f file. I don't want to run the whole pipeline (including alignment). I was trying to manually generate a hd5f file to input to cellranger reanalyze. Is this possible?

The direct answer is yes. If you have never worked with hdf5 you can start  here.

But, cellranger outputs several hdf5 files. You can extract the count matrix from the  "/outs/filtered_gene_bc_matrices_h5.h5" file.

They provide a lot of information on how to use the hdf5 format so their their documentation is a good place to start.

In addition: you can also load the data into Seurat, monocle or scanpy and follow their tutorials to analyze the data.


How to remove sequences from a fasta file using a sequence ID list which contains a space within the id? I have a fasta file that contains sequence reads and sequence id file that needed to be removed from the fasta file. I have done this earlier, but since id contains a space my piece of code is not working and I don't know how to change it since I'm still learning. fasta file looks like this.
>m64041_200717_231916/74/12723_24868 id=30
CCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCA
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/105/20691_65844 id=488
AGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCAAGACTTTAT
>m64041_200717_231916/108/17414_67048 id=4956
TGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

IDList file is like this.
m64041_200717_231916/74/12723_24868 id=30
m64041_200717_231916/108/17414_67048 id=4956
m64041_200717_231916/105/20691_65844 id=488

I need to remove these sequences from that fasta file. I have a code I used earlier. But it is not working since id contains a space (I think). Previous code is like this.
# use as follows
# ./filter_seq.bash <fastafile> <list_of_ids_you_want_to_remove>
# ./filter_seq.bash all.fasta bacterial_IDlist

grep ">" $1 | sed 's/>//g' | sort | uniq > all_del
sort $2 | uniq > query_del
comm -23 all_del query_del | sort | uniq > filtered_del
perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' filtered_del $1 > filtered.fasta
rm *_del

I would be grateful if someone can help me.

Install (pip) biopython if u don't have it.
from Bio import SeqIO
with open('/path/ids.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('/path/output.fa', 'w') as fout:
    for record in SeqIO.parse('/path/sequences.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')


How is consensus alignment for OLC assembly usually implemented? Background

In C++ I want to implement a simple class that performs overlap layout consensus assembly but I can not figure out the most logical data structure to use for the consensus alignment step.

This is the step where all of the reads from the overlap graph are aligned and the consensus is called (slide 48). An example alignment is below. Asterisks above the column indicate some disagreement among all reads. Alignment position is shown above and below for convenience.

                         *    * *         *
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012
read 1   TAGATTACACAGATTACTGAT TGA
read 2      ATTACACAGATTACTGACTTGATG CGTAAACTA
read 3                   tTGACTTcATGGCGTAA CTA
read 4                     GACTTGATGGCGTAA CTAAG
read 5                           ATGGCGTAA CTAAGATCG
     pos           111111111122222222223333333333444
         0123456789012345678901234567890123456789012

Question:

Q1: What is the most logical way to implement the consensus alignment step in a OLC assembler?

Q2: What are the pros and cons of this and other possibilities?

Ideas

Idea - Vector of vectors. One column per position.

My initial idea for a data structure to store the alignment is to use one vector per read to store the read's character at that position. We would end up with j arrays (number of reads) all of length i (the total length of the overlapped alignments). This, however, would be space-inefficient if the locus were large and many positions in the vectors would contain a null value. In the example alignment above the first 23 positions of read 5 would be null values.

Idea - Graph. One node per position.

In another implementation idea, each position in the alignment would be represented by a node in a directed graph. Each node would have a property keeping track of the nucleotide counts at each position and might keep the read id's stored as well if there is a good reason to do so.

For short reads, mainstream assemblers don't have a consensus step. The contig path directly spells the consensus sequence.

For long reads, the most popular way to derive a consensus is via partial-order alignment (POA). This is used by error correctors (similar to consensus) like pbdagcon and falcon_sense, and consensus tools like canu, racon, wtdbg2 and nanopolish. Basically, you start with an empty directed acyclic graph (DAG), align each read in turn to the graph and augment the DAG by inserting unaligned segments. When you align all reads onto the DAG, the most weighted path is taken as the consensus. You should read Jared Simpson's blog on this topic. It is much more detailed than "basically". You should also read the papers on all the mentioned tools above if you want to implement by yourself.

An alternative approach is Gene Myers' realigner (PDF here). It is not often used, though.


Retrieving NCBI Taxa IDs from refseq or GenBank assembly accession I have about 10,000 genome files all named by either refseq or genbank accession number, do you know if it's possible to convert these numbers to the corresponding NCBI taxon ID or species?

for example:

GCA_000005845.2 to 79781

In the case of E.coli

Edit:

The file names look like this:

GCF_000337855.1_ASM33785v1_protein.faa.gz
GCF_000091125.1_ASM9112v1_protein.faa.gz
GCF_000184535.1_ASM18453v1_protein.faa.gz

My operating system is Ubuntu 16.04

Turns out I'd already written some code that did it, human memory is a funny thing.

This uses biopython to split the field description to where the species is. May not work for all NCBI files, but seems to work on most.

import Bio
from Bio import SeqIO
from Bio import AlignIO

for record in SeqIO.parse (FILE, "fasta"):
    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]


Determining position of side chain hydrogen in glycine residues from coordinates of backbone atoms I have a dataset of pdb-format protein structure decoys/models generated by tertiary structure prediction programs. Many of these are missing hydrogen atoms.

Specifically, I'd like to determine the rough location of side-chain hydrogens in glycine residues, given the CA, C, N, and O coordinates. Is this possible? I don't need too much precision - I just need the general direction relative to the rest of the residue. More specifically, I want a vector originating at the CA that points roughly in the direction of the side chain.

I've looked at software for generating hydrogen atoms from pdb structures, however, I need to know which ones are glycine side chains, and was wondering if there might be a simpler way to predict hydrogen locations given coordinates of backbone atoms.

Glycine has a single hydrogen atom as its side chain:

All the six bond angles with the CA atom in the middle are about 109°
(C-CA-N, C-CA-HA3, C-CA-HA2, N-CA-HA3, N-CA-HA2 and HA3-CA-HA2 using the CCD naming convention). This defines the rough direction of the hydrogen atoms.

But which of the two hydrogens is the side chain?

If it is a left-handed protein, I'd assume that it's the hydrogen pointing toward the screen, similarly to the side chain in L-alanine:

The images are taken from WikiMedia Commons:
glycine
and L-alanine.


Question regarding the function pairwiseAlignment in R I'm a bit confused and would appreciate your assistance, please.

I'm trying to use the pairwiseAlignment function from the Biostrings package to align two protein sequences.

I'm not quite sure where to put my reference sequence, the sequence I want to compare to. Is it under "subject" or "pattern"?

Also, is the score I get the edit distance? Is it the alignment score? Or is it the the alignment length? When I use local alignment, will the score be higher the more the sequences resemble each other? And when I use global alignment, will the score be lower the more the sequences resemble each other?

According to the help page, it seems more logic to put the reference sequence in subject, and your query sequence in pattern, but it won't really matter for the outcome (except that the pattern sequence will be on top of the alignment).

For your second question:

  the scores for the optimal pairwise alignments is returned


Installing Rostlab's PredictProtein I am trying to install Rostlab's PredictProtein in a Virtual Machine as detailed here https://rostlab.org/owiki/index.php/How_to_Run_The_PredictProtein_Machine_Image_with_VirtualBox:

  Download PredictProtein Machine Image to your local drive. You will need the Bio-Linux-based 'VMWARE/VMDK' format version.

  Choose the Debian-based version in case you do not need a graphical
  user interface. This image is only half the size of the
  Bio-Linux-based image.

  Download a database image to your local drive. We recommend the 'raw' format version.
  Create New Virtual Machine: Machine -> New
  Name your virtual machine as 'PredictProtein'
  Choose Linux and Debian as operation system
  Locate and select the VMware .vmdk file as your boot disk
  Give at least 2GB of RAM to your virtual machine
  Click 'Finish'
  Boot your new virtual machine.
  If the message "The disk drive for /mnt/local-storage is not ready yet or not present" appears, press 's'.

The problem comes at this step:

  Open a terminal and enter:

sudo mkfs.ext4 -F -m 0 -L local-storage /dev/sdb
sudo mount -a
sudo mkdir /mnt/local-storage/rostlab-data && cd /mnt/local-storage/rostlab-data
sudo mount -a
sudo mount /mnt/rostlab-data
sudo tar -xvJf /mnt/rostlab-data/rostlab-data.txz
sudo umount /mnt/rostlab-data && cd

When I run sudo mount /mnt/rostlab-data, I get:

mount: special device LABEL=rostlab-data does not exist

How should I proceed?

Additional notes -

The output for

cat /etc/fstab

is

# /etc/fstab: static file system information.
#
# Use 'blkid -o value -s UUID' to print the universally unique identifier
# for a device; this may be used with UUID= as a more robust way to name
# devices that works even if disks are added and removed. See fstab(5).
#
# <file system> <mount point>   <type>  <options>       <dump>  <pass>
proc            /proc           proc    nodev,noexec,nosuid 0   0
# / was on /dev/sda1 during installation
UUID=9a9418ff-48f9-421a-b25a-82161839351d / ext4 errors=remount-ro 0 1
# swap was on /dev/sda5 during installation
UUID=a96757f4-8cc7-410f-95ee-b5ae642eb9ed none swap sw  0       0

LABEL=local-storage /mnt/local-storage auto defaults,nofail 0   2
LABEL=rostlab-data /mnt/rostlab-data auto defaults,noauto,ro 0  0

/mnt/local-storage/rostlab-data /usr/share/rostlab-data none bind,nofail 0 0

The output for

blkid

is

/dev/ramzswap0: TYPE="swap"
/dev/sda1: UUID="9a9418ff-48f9-421a-b25a-82161839351d" TYPE="ext4"
/dev/sda5: UUID="a96757f4-8cc7-410f-95ee-b5ae642eb9ed" TYPE="swap"
/dev/sdb: LABEL="local-storage" UUID="df841028-8dd1-4cd2-b1bb-6ae1f6887510" TYPE="ext4"

It looks like the image you downloaded is missing the partition with the label rostlab-data. Your blkid output shows that you have the local-storage partition, but not the rostlab-data which is automatically added to your /etc/fstab.

I am afraid you'll have to contact the developers and let them know. Sounds like a bug.


Stranded vs. unstranded library preparation protocols in RNAseq I've been reading this paper lately:

Sailfish Enables Alignment-Free Isoform Quantification from RNA-Seq Reads Using Lightweight Algorithms

I don't really understand the second paragraph under Quantification in the Online Methods section. I quote the author:

  In an unstranded protocol, sequencing reads, and hence the k-mers they contain, may originate from transcripts in either the forward or reverse direction.

As far as I understand we are sequencing cDNA here and it is a single stranded molecule so all the reads will have the same direction (unlike the DNA where the reverse complement strand is also there), so I don't see the author's point. Can someone clarify please? Thanks.

Short answer: cDNA is double-stranded.

I'll try to explain the concepts surrounding this with a nanopore read, where the cDNA component of the read (from the mouse genome) is indicated in a light grey, and the adapter components are either coloured or dark grey:

The cDNA library from this RNA was prepared by adding a polyT adapter (binding to the polyA sequence of the RNA), and carrying out reverse transcription to the start of the transcript. A strand-switching process was then used to prime and extend along the opposite strand of the DNA, creating a double-stranded full-length cDNA template. This template was then amplified using PCR to increase cDNA concentration and add in barcoding sequences.

You may notice that this particular read has a polyT head, rather than a polyA tail, attached to the VNP primer in the forward orientation. I can tell by this feature that this read comes from the reverse-complement strand of a DNA transcript, and so I know that the actual RNA transcript sequence is the reverse-complement of this read. Both the forward and reverse-complement sequences could potentially come through in the sequencer, because the adapters for sequencing are attached to both sides of the cDNA template, but the transcript direction can still be determined by looking at the orientation of the adapter sequences. When mapping this sequence to the genome, only a single strand (and its associated set of kmers) needs to be considered.

When considering what happens after mapping to the genome, there are four possible match situations:

Adapter is in reverse orientation, genome is in reverse orientation
Adapter is in forward orientation, genome is in reverse orientation
Adapter is in reverse orientation, genome is in forward orientation
Adapter is in forward orientation, genome is in forward orientation

In the case of the library preparation above, situations 2 and 3 correspond to a transcript that is encoded on the forward strand of the genome, whereas situations 1 and 4 correspond to a reverse-encoded transcript. Different library preparation methods could have a different way of encoding transcript direction. See here for more details on how these match situations can be used to produce strand-specific BAM files.

If I were to use an unstranded protocol (e.g. generating a cDNA library through other means, and attaching adapters randomly to the forward or reverse transcript strand), then the adapter direction wouldn't give me information about which strand this sequence came from. I would need to consider both strands when mapping (or more correctly, when identifying the native direction of the originating transcript), because the adapter sequence has been decoupled from the transcript direction.

On Illumina systems, this strand direction calculation is performed via the chemistry on the flow cell, and the adapters are [usually] not sequenced by the sequencer at all. Adapters from one direction bind to complementary sequences on the flow cell and are extended to produce the first read, and adapters from the other direction bind to produce the second read. Even with a stranded protocol, adapters from both directions are bound (because the cDNA template is double-stranded), but the binding order (as first or second read) can be used to determine the native direction of the originating transcript.


BinaryCIF vs MMTF formats, which one to choose? Both file formats are advertised as more memory/parsing efficent than old PDB format and its successor mmCIF format, but I can't understand the context (I'm a developer, but pretty new to the field of bioinformatics).
Is BinaryCIF a later iteration of MMTF concepts/ideas (so it's a new competing format) or is the purpose of the two formats different ?
Refs :

BinaryCIF: https://github.com/dsehnal/BinaryCIF
MMTF: http://mmtf.rcsb.org/index.html

So binaryCIF is newer and is heavily inspired by MMTF.
NGL, the most "loved" JS library (IMO) to show protein fetches proteins as MMTF, while Mol* uses bCIF —Mol* has the official repo of the bCIF format. Alex Rose, developed the former at the RCSB PDB, and the latter in collaboration with the PDBe. He also worked with Antony Bradley, the main author on the MMTF format, who left the RCSB PDB, some 4+ years ago. So there is definitely no competition there.

The MMTF download of the RCSB works well: https://mmtf.rcsb.org/download.html
The MMTF format did not catch on —okay, mmCIF is only slowly replacing PDB, which is still the only accepted format for analyses—, hence the need for a bCIF format
The libraries to do the conversion aren't without glitches
In the distant future bCIF will be the sole binary format
mmCIF <--> bCIF conversion is easier

However, I really ought to stress that mmCIF is a better choice to bCIF or mmTF. And even then many tools still hate mmCIFs and use PDB format.
If memory footprint is not really a problem, then stick with mmCIF. If you are running a webserver, you may be sending data as gzipped making this a moot point.


Adjusting phenotypes by regressing out covariates I'm trying to use the bfGWAS tool, which analyses GWAS data and integrates functional annotations to identify casual SNPs (paper and
github). In the user manual, it states:

  We recommend first regressing covariants out from the original
  phenotypes, and then provide bfGWAS the corrected phenotypes (i.e.,
  residuals from the regression model with covariates).

It is unclear to me what this means and how to go about determining the corrected phenotypes. Typically, in an association test with covariates, a regression would be performed using a model like

pheno ~ snp + covariate

I could then get the residuals from this model. However, I'm not certain how to go about it for all SNPs, and then adjust the phenotypes. I have case-control data, so the phenotypes are binary.

Instead of the "observed phenotype", it's asking for a "corrected phenotype", which is the residuals from pheno ~ covariates. In R, one would get that as follows:

fit = lm(phenotype ~ covariates, data)
corrected_phenotype = resid(fit)

Of course covariates is just a place hold for something like age + gender or something along those lines.


Annotatation of transposable elements / transposons / SINE sequences in human I need to annotate human genome with different types of information, amongst them: transposable elements.
I read about repeatmasker, but what I understood is that it estimates the probability of a sequence to be a transposon based on the sequence and not based on curated database.
My university does not have access to RepBase (GIRI). I was hoping finding some possibility with biomart but it seems not possible.
What would be a way to access curated annotation of all kinds of transposable elements (SINE, Alu seq..) within a given genome/genome coordinates?
My goal is to estimate the link between some measurements we have made on the genome (binned) ( and Shannon entropy is one of them, and probably the less interesting regarding TE). Eventually, I'd like to be able to estimate that measurement X has a correlation with the number of Alu repeats present in the neighbourhood.
What would be the best way to achieve this?

repeatmasker is meant to annotate/mask a genome given a library of known transposable elements. From the project description, I assume you have not done any sequencing on your own. Then, instead of annotating, I think you should find already computed annotations of TEs.
The alternative to subscription-based repbase is dfam. And Just looking at their webpage they seem to have an annotation of all individual TEs in hg38. I bet my shoes there will be plenty of TE annotations in human in genome databases too.


Integrative analysis of omics studies using machine learning I would like to use public omics datasets (ChIP-seq, RNA-seq, and ATAC-seq) from different studies to do an integrative analysis as follow:

Normalise samples, within each type of omics, from different public datasets.
Convert the normalised values into a uniform scale to make the comparison between ChIP-seq, RNA-seq and ATAC-seq possible.
Feed the normalised uniformed values into machine learning to infer one feature (e.g. RNA expression) from other features (e.g. TF or histone marks ChIP-seq).

How could I integrate these data-sets to predict one group of variables from another one?

The steps you describe are correct. For step 2 it is usually normalized to mean 0 and variance 1. However the "machine learning" part is important.

Having several samples being technical replicates will make the integration task easier. However, you have too few samples to make any good prediction. At most I would describe it as an exploratory analysis.

I recommend using regularized canonical correlation analysis for this kind of multi-omics analysis. The method seeks which are the features of one dataset that correlate (or covariate, depending on a parameter) with the variables of other datasets. Then you could use those features to predict the position of the other.

There are several implementations of this methods, the original in the RGCCA package, and mixOmics for other procedures with cross-validations also STATegRa has other related methods.


Using ONT MinION, why is that we cannot get a full length DNA read? Here, https://www.youtube.com/watch?v=E9-Rm5AoZGw, ONT claims they can read non-fragmented DNA and that with the hairpin like structure the double strands can be read at once. Then why is that MinION reads are of variable in length and not the length of a whole DNA? BUt here:https://www.youtube.com/watch?v=e9wtAIvPPxY, ONT claims that the read lenght can be configured, what does this mean?

The MinION should be able to sequence DNA fragments entirely, regardless of the size, if you can get them to the pore intact. That is the problem: during DNA extraction and library preparation you will introduce breakage.


Duplicate long hits from PSI-BLAST I had a protein Refseq ID and I PSI-BLASTed this sequence against Refseq database. We all know that the Refseq is a Reference sequence database and it shouldn't have redundancy. After BLASTing my sequence, at first iteration I got 1000 hits and among them there were a lot of redundant sequences!

My sequence had 241 amino acids and I found a lot of sequences with 100% identity, 100% cover and 0 E-value exactly the same as my sequence but with different IDs. All of these IDs were from RefSeq! In other iterations and after adjusting format options, I got this redundancy with other sequences from other species. My sequence is related to a chain of a multichain protein (E.coli fumarate reductase)

For example, when using NP_418578 as a query, I found WP_078165098.1, WP_064226696.1, WP_062863447.1, WP_001401474.1 and other that were identical.

I want to know what is wrong with Refseq. Is it really a Reference sequence database? If it is, where does this redundancy come from and why?

This is what is known as a feature, not a bug. Note that your identical proteins all have accessions starting with WP_. These are special "non-redundant proteins". Many sequences—particularly bacterial sequences—are identical between various different species so having a separate RefSeq entry for each of them would be inefficient. Therefore, RefSeq combines multiple such proteins into a single WP_ record. This is documented here (emphasis mine):

A new type of RefSeq protein record which represents non-redundant protein sequences was introduced in mid-2013. This record type was introduced to address a growing issue with redundancy in the Prokaryotic RefSeq protein dataset that coincided with a significant increase in bacterial genome submissions from individual isolates and closely related bacterial strains. For example, a large number of high-quality bacterial genomes may be submitted during a disease outbreak. The submitted sequences may reflect pathogen evolution during the course of the outbreak but the majority of the encoded proteins from these genomes may be identical to each other. As RefSeq includes these genomes, per community requests, this resulted in increased redundancy. By representing identical proteins using a single non-redundant protein accession number (with the prefix 'WP_'), redundancy in the database is significantly reduced.
[ . . . ]
Because a non-redundant protein sequence may be found in RefSeq genomes from multiple species, the organism information provided on the protein record reflects the lowest-common taxonomic node ranging from the genus species level to super-kingdom. A non-redundant protein record that provides organism information at the level of a genus, family, or even super-kingdom does not mean that the protein is found in all RefSeq genomes below that taxonomic classification. It only indicates that the protein is found in more than one genome of different species for which the genus, family, or super-kingdom classification is the lowest common taxonomic node.

So, your query was NP_418578.1, anaerobic fumarate reductase catalytic and NAD/flavoprotein subunit from E. coli strain K-12, substrain MG1655. The first thing to notice is how specific that is. This is the protein found from one specific substrain of one specific strain of one specific bacterial species. It is reasonable to expect that there will be identical sequences from many, many closely related species. Both from, most probably, all other strains and substrains of E. coli and from other, related bacteria.
Now, the specific sequences you mention are actually slightly different and not 100% identical. Below is a multiple alignment of NP_418578.1 and the 4 WP_ sequences you mentioned. Note that each of the 5 entries is slightly different. Each has one residue that differs from the rest. Look for the : in the identity line, there are 4 : and all others are * (I am only showing the relevant alignment blocks here and have removed those where all 4 sequences were identical):
WP_001401474.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_062863447.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_064226696.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
WP_078165098.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAIAQDH
                    *******************************************************:****

[ . . . ]

WP_001401474.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_062863447.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_064226696.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHIRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
WP_078165098.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ***********************************************:************

[ . . . ]

WP_001401474.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_062863447.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKISQAFWHEWR
WP_064226696.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
WP_078165098.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    **************************************************:*********

WP_001401474.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGIDPVKEPIPVRPTAHYTMG
WP_062863447.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_064226696.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
WP_078165098.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    *****************************************:******************

[ . . . ]

Your sequence (NP_418578.1) is only identical to one WP_* multi-species sequence, WP_001192973:
WP_001192973.1      MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
NP_418578.1         MQTFQADLAIVGAGGAGLRAAIAAAQANPNAKIALISKVYPMRSHTVAAEGGSAAVAQDH
                    ************************************************************

WP_001192973.1      DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
NP_418578.1         DSFEYHFHDTVAGGDWLCEQDVVDYFVHHCPTEMTQLELWGCPWSRRPDGSVNVRRFGGM
                    ************************************************************

WP_001192973.1      KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
NP_418578.1         KIERTWFAADKTGFHMLHTLFQTSLQFPQIQRFDEHFVLDILVDDGHVRGLVAMNMMEGT
                    ************************************************************

WP_001192973.1      LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
NP_418578.1         LVQIRANAVVMATGGAGRVYRYNTNGGIVTGDGMGMALSHGVPLRDMEFVQYHPTGLPGS
                    ************************************************************

WP_001192973.1      GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
NP_418578.1         GILMTEGCRGEGGILVNKNGYRYLQDYGMGPETPLGEPKNKYMELGPRDKVSQAFWHEWR
                    ************************************************************

WP_001192973.1      KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
NP_418578.1         KGNTISTPRGDVVYLDLRHLGEKKLHERLPFICELAKAYVGVDPVKEPIPVRPTAHYTMG
                    ************************************************************

WP_001192973.1      GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
NP_418578.1         GIETDQNCETRIKGLFAVGECSSVGLHGANRLGSNSLAELVVFGRLAGEQATERAATAGN
                    ************************************************************

WP_001192973.1      GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
NP_418578.1         GNEAAIEAQAAGVEQRLKDLVNQDGGENWAKIRDEMGLAMEEGCGIYRTPELMQKTIDKL
                    ************************************************************

WP_001192973.1      AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
NP_418578.1         AELQERFKRVRITDTSSVFNTDLLYTIELGHGLNVAECMAHSAMARKESRGAHQRLDEGC
                    ************************************************************

WP_001192973.1      TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
NP_418578.1         TERDDVNFLKHTLAFRDADGTTRLEYSDVKITTLPPAKRVYGGEADAADKAEAANKKEKA
                    ************************************************************

WP_001192973.1      NG
NP_418578.1         NG
                    **

So, in summary, RefSeq will combine multiple identical sequences into a single WP_* multi-species accession. You should therefore expect to find one 100% identical WP_* sequence for your query and multiple, almost identical WP_* entries. And that's precisely what you see here.


Download multiple SRA files I want to download all SRA file from the following project. Is there a method to download all the SRA files at the same time?

A quick look at your link tells me the SRR numbers run from SRR837819 to SRR837856. You can use fastq-dump from the sratoolkit, and make a for loop around it in bash.

Something like this should work:

for (( i = 19; i <= 56; i++ ))
  do
  fastq-dump --accession SRR8378$i
done

After reading Devon Ryan's answer, I realize that you asked for SRA files instead of fastq. This can also be done with prefetch of the sratoolkit:

for (( i = 19; i <= 56; i++ ))
  do
  prefetch SRR8378$i
done


Python: How to write duplicate sequences removed from fasta file to new file I currently am using this code to remove duplicate sequences from the fasta file. However, I would also like to write a new file with only the removed duplicates as well as a count for how many times they appear. Please let me know if you have any suggestions, I am hoping to keep it written in one py file, but I am not opposed to running a second file to do so. I attended to write the appended strings to a new file but it looked slightly off...
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   if str(record.seq) not in seen:
       seen.append(str(record.seq))
   records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")
```

Please do not reinvent the wheel - seqkit rmdup does what you want.
From the website:

rmdup
Usage
remove duplicated sequences by id/name/sequence

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case

How does a -D file look like, you ask:
$ zcat hairpin.fa.gz \
    | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2


Error compiling HMMER v3.2.1 I tried to compile HMMER v3.2.1 in Cygwin using

./configure

followed by

make

or

make install

but failed. I got this config.log file to troubleshoot and I suspect the solution is in this line:

configure:2908: error: no acceptable C compiler found in $PATH

What does it mean? Is it really all I need to address? I followed the installation steps in this document.

I would strongly recommend to use Bioconda instead of compiling. It is one of the channels ("distributors of packages") in Anaconda Cloud which is basically a package source for conda.

In short - conda allows you to create virtual environments isolated from system and to install various precompiled software packages with one command (so, by the way, you can easily use multiple software versions on one system). Also, Bioconda is just one channel aiming on various bioinformatics software, but many more packages can be found in aforementioned Anaconda Cloud.

Tutorial how to start using Bioconda is on its homepage, and it is really easy (far more easier than to struggle with compiling).

Bioconda HMMER package

P.S. Of course there are cases in which you need a specific build of given software. Than it is really needed to compile manually.


Generating basic dna sequences in R Looking for some guidance.

In an introductory bioinformatics class I am tasked with doing the following in R:

Create the sequence ACTGACTG
Generate a random string of nucleotides with equal probabilities.

Are there certain functions that I should be looking at? No library or package or set of functions was provided so I am completely lost. I am assuming there must be functions in R that do these things that I need to manipulate, but I cannot seem to find any resources explaining that. I looked at seq() (thinking it might generate DNA sequences), but it only seems to generate sequences of integers. I also looked at randDNA(), but it is random, how could I make it generate the exact string in question 1).

Any direction will be greatly appreciated.
Thanks.

Without giving away too much for your homework check out the sample() function in R. Think about how you can pass DNA bases to this function to generate random sequences.


How to simulate "base error rate" in art_illumina? I'd like to simulate 10% sequencing error using art_illumina. The simulator doesn't have a parameter that I can just give the 10%, but it has this:

-qs  --qShift   the amount to shift every first-read quality score by
-qs2 --qShift2  the amount to shift every second-read quality score by
                NOTE: For -qs/-qs2 option, a positive number will shift up
                quality scores (the max is 93) that reduce substitution sequencing
                errors and a negative number will shift down quality scores that
                increase sequencing errors. If shifting scores by x, the error
                rate will be 1/(10^(x/10)) of the default profile.

Heng Li's wgsim has an option for base error rate. Can I do the same for art_illumina?

Q: Does the "error rate" in the description mean sequencing error? If I set both qShift and qShift2 to 10, 1/(10^(x/10)) == 0.1, where x is 10, does that mean 10 is the value I should give to the simulator?

You need to specify an error profile to do this. However, you can shift quality scores of an existing profiles to get 10% error, which is easier than defining a complete profile.

For example, Illumina HiSeq 1000 has the error rate of 0.1% and you want 100 time more errors, so you should shift down the quality scores by 20 by setting both qShift and qShift2 to -20. The following code will simulate paired end reads with 10% error.

art_illumina -ss HS10 -i reference.fa -qs -20 -qs2 -20 -sam -p -l 100 -f 20 -m 200 -s 10 -o output


How to convert ChiP-seq track in .bam file format into an array at fixed bin resolution? Is there a out of the box to convert a .bam file to a .wig file at specific resolution/bin size?

I was thinking to use igvtools, but I am not sure how it works:

igvtools count -z 5 -w 25 -e 250 alignments.bam  alignments.cov.tdf  hg19

What is the diffrence between -w,-z,-e?

-z num  Specifies the maximum zoom level to precompute. The default
               value is 7 and is sufficient for most files. To reduce file
               size at the expense of IGV performance this value can be
               reduced.

 -e, --extFactor num

The read or feature is extended by the specified distance in bp prior to counting. This option is useful for chip-seq and rna-seq applications. The value is generally set to the average fragment length of the library minus the average read length.

-w, --windowSize num

The window size over which coverage is averaged. Defaults to 25 bp.

Using deepTools to get a bin size of 10 and 20 threads:

bamCoverage -of bedgraph -bs 10 -b input.bam -p 20

You'll probably prefer the default bigWig format if you need random access or to use something like IGV.

Regarding the options to igvtools, -w is equivalent to -bs/--binSize in deepTools and denotes the resolution of the output data. -e (it's the same option in bamCoverage) denotes whether and how much to extend a read in its 3' direction. This done such that the output represents the sequences fragments. -z is specific to things like bigWig files (and apparently tdf files), where there are precomputed values that IGV will load as you zoom out (the various resolutions are "zoom levels", which is parlance taken from the bigWig format).


samtools depth print out all positions I am trying to use samtools depth (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage at every position:

cat GRCh38.karyo.bed | awk '{print $3}' | datamash sum 1
3088286401

I would like to know how to run samtools depth so that it produces 3,088,286,401 entries when run against a GRCh38 bam file:

samtools depth -b $bedfile -a $inputfile

I tried it for a few bam files that were aligned the same way, and I get differing number of entries:

3087003274
3087005666
3087007158
3087009435
3087009439
3087009621
3087009818
3087010065
3087010408
3087010477
3087010481
3087012115
3087013147
3087013186
3087013500
3087149616

Is there a special flag in samtools depth so that it reports all entries from the bed file?

If samtools depth is not the best tool for this, what would be the equivalent with sambamba depth base?

sambamba depth base --min-coverage=0 --regions $bedfile $inputfile

Any other options?

You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.

bedtools genomecov -d -ibam $inputfile > "${inputfile}.genomecov"

You can also provide a BED file if you just want to calculate in the target  region.


What Ensembl genome version should I use for alignments? (e.g. toplevel.fa vs. primary_assembly.fa) When you look at all the genome files available from Ensembl. You are presented with a bunch of options. Which one is the best to use/download?

You have a combination of choices.

First part options:

dna_sm - Repeats soft-masked (converts repeat nucleotides to lowercase)
dna_rm - Repeats masked (converts repeats to to N's)
dna - No masking

Second part options:

.toplevel - Includes haplotype information (not sure how aligners deal with this)
.primary_assembly - Single reference base per position

Right now I usually use a non-masked primary assembly for analysis, so in the case of humans:
Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

Does this make sense for standard RNA-Seq, ChIP-Seq, ATAC-Seq, CLIP-Seq, scRNA-Seq, etc... ?

In what cases would I prefer other genomes? Which tools/aligners take into account softmasked repeat regions?

There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats are (we never make use of this for our *-seq experiments, but it's there in case we ever want to).

For primary vs. toplevel, very few aligners can properly handle additional haplotypes. If you happen to be using BWA, then the toplevel assembly would benefit you. For STAR/hisat2/bowtie2/BBmap/etc. the haplotypes will just cause you problems due to increasing multimapper rates incorrectly. Note that none of these actually use soft-masking.


What is the difference between SAM mapping quality and Blast E-value? Blast reports E-values, but short-read mappers report mapping qualities. Are they the same thing? Can they be converted to each other? If not, why blast doesn't report mapping quality while short-read mappers do not report E-values?

The E-value and the mapping qualities are two very different things.

The E-value is "a parameter that describes the number of hits one can 'expect' to see by chance when searching a database of a particular size". More details can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect

The mapping quality is an attempt to estimate the probability that a given base from a given read is mapped correctly to a particular place in the reference genome. Different aligners calculate mapping qualities very differently, so there is no simple way to compare them across aligners. But in general the mapping quality will consider things like: (i) the quality of the base call; (ii) the repeat structure of the reference; (iii) the alignment algorithm; (iv) whether the read has a mapped pair; (v) anything else the author of the mapping software thought might help.

E-values and mapping qualities are, therefore, measuring two fundamentally different properties with very different uses.


Find specific functional genes from a list? I have compiled a list of genes that are distinct from 2 different conditions.
From this list (of about 600 human genes), I want to find out if any of the genes in this list have acetyl transferase function i.e. are acetylases.

Is there an online platform where I could say put in this list of genes and the output I can select genes which have this property/function?

I.e some sort of gene annotation for functions for a list of human genes?

BioMart. Filter by your list of genes and a GO term linked to acetyl transferase activity.


How to choose docking score cutoff? I made a virtual screening of a large database (~5M) and want to filter ligand structures for further more accurate screening (from Schrodinger HTVS to SP resolution). My next step is to define a cutoff docking score, but it seems arbitrary to choose one. Are there any established methods for filtering the ligands out?

It is utterly arbitrary and depends on your resources.
We can all agree your cutoff would be a lot lower during the ordering step than at followup in silico analyses, but even then if you were ordering the top N at \$100-250 each, should you spend \$1,000 or \$10,000?
In the in silico step you have real time and CPU-time to factor in making this call even blurrier.
But there are some things to put those scores into context:

The units are kcal/mol, i.e. these are predicted ∆∆G of binding. Here is an explanation of what does a kcal/mol mean for docking — 1 kcal/mol is a hydrogen bond or the average collision energy of water at 37°C.
If you have a crystal structure with a known ligand, re-docking and cross-docking that would give you a benchmark of what you should expect ——keeping my mind that natural ligand are horrid and don't obey Lipinski's rule. The natural ligand will wander and to assess how good the docking went RMSD is used, but for this I mean looking at its GlideScore as if it were any other compound.
Rescore with a different tool —Gold etc.

Ligand efficiency
An important thing to remember is that ∆∆G is size dependent —a huge peptidomimetic will score more than a smaller molecule that hits the hotspots perfectly. Therefore, ligand efficiency (LE), ∆∆G over heavy atoms (reviewed in detail here) is a good metric to prevent huge molecules dominating spuriously. If your library has a variety of different sizes, it is best to sort by LE than ∆∆G alone. In fact, if you get a few small hits you can merge/link them or expand on them, whereas if you get nothing because you gambled on larger hits you have to start from square one. I would argue that this is true for repurposing screens and the low success of repurposing VLS despite the enthusiasm is due to the neglect of this.


Marking optical or PCR duplicates with picard vs. samtools flagstat I am trying to identify the best library prep method for noninvasive prenatal test samples, to be sequenced with Illumina Novaseq.

One metric that I am evaluating is the number/ percentage of PCR or optical duplicates. In order to achieve that I have done the following steps:

aligned the fastq reads against the reference genome with bwa mem
the .sam output from bwa mem was converted to bam, sorted and indexed
the duplicates were marked with the command MarkDuplicates from picard

Then if I call samtools flagstat on the sorted bam file which had the duplicates marked with picard I get:

26595942 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
466 + 0 supplementary
1636809 + 0 duplicates
24969064 + 0 mapped (93.88% : N/A)
26595476 + 0 paired in sequencing
13297738 + 0 read1
13297738 + 0 read2
21102678 + 0 properly paired (79.35% : N/A)
24049752 + 0 with itself and mate mapped
918846 + 0 singletons (3.45% : N/A)
462212 + 0 with mate mapped to a different chr
182138 + 0 with mate mapped to a different chr (mapQ>=5)

Which sounds a bit strange to me, as no duplicates were found.

However, the MarkDuplicates command also outputted a metrics file (dups/dupsMetrics.txt) that did identify some duplicates. This was achieved with the following command:

java -jar picard.jar MarkDuplicates I=alignment/sample.sorted.bam O=alignment/sample_markDup.bam M=dups/dupsMetrics.txt

I see that the value under the column PERCENT_DUPLICATION in dups/dupsMetrics.txt is 0.065555 (given the column name, I would expect that the value is indeed 0.06%, and not 6.55%)

What would be the explanation between the discrepancy in the output of samtools flagstat and MarkDuplicates. Is it reasonable to expect that the value displayed by samtools flagstat is rounded to zero, given the low percentage which was observed by MarkDuplicates (0.06%)?

There are duplicates, in this line:

1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356

According to samtools documentation for flagstat:

  Provides counts for each of 13 categories based primarily on bit flags
  in the FLAG field. Each category in the output is broken down into QC
  pass and QC fail. In the default output format, these are presented as
  "#PASS + #FAIL" followed by a description of the category.

So you can usually take whatever is in the first column


Getting all NNI trees of a parsimony tree I have an aligned protein sequence file which I have been using for reconstructing a parsimonious tree. I am currently using NNITreeSearcher._get_neighbors method from Biopython 1.72 but it's way to find one best scored tree only.

def _nni(self, starting_tree, alignment):
    """Search for the best parsimony tree using the NNI algorithm (PRIVATE)."""
    best_tree = starting_tree
    while True:
        best_score = self.scorer.get_score(best_tree, alignment)
        temp = best_score
        for t in self._get_neighbors(best_tree):
            score = self.scorer.get_score(t, alignment)
            if score < best_score:
                best_score = score
                best_tree = t
        # stop if no smaller score exist
        if best_score >= temp:
            break
    return best_tree

I need to generate all NNI-neighbours of the trees (two trees are NNI-neighbours if one can be transformed into another by one nearest neighbour interchange operation).

I want to know how can I get total number of trees generated before the best tree is reached.
Is it possible to return all the trees generated before the best tree is reached, so that the all/required number of trees(other than the best tree) may be printed?

As I've studies the 'TreeConstruction' module of Biopython, I came to know that '_get_neighbors' method is the one that is used to get all neighbor trees of the given tree. This is all I've concluded from my findings. If anyone wants to add more, that'll be positively appreciated.


Download or create VCF from Human Genome Diversity Project (HGDP) I can't seem to find HGDP data in VCF format.
Does anyone know where I can download it from?
Perhaps there a Snakemake pipeline for calling VCFs from this data?

I presume you are looking for the recently sequenced WGS data rather than the genotype data? If so, the different files can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/
Each chromsome has a seperate vcf file with a corresponding .tbi index file.
The gvcfs can be found here: ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/gVCFs/
These files are more compressed and smaller in size than the regular vcfs.


what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets? What is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets nowadays?

About 5 years or so ago, it used to be the latest release of the CSHL HapMap 12 populations that were part of the 1000 genomes project. These were:

CHB CHD MEX GIH TSI LWK ASW MKK HCB JPT CEU YRI

For example, the EnsEMBL project currently has these 12 populations available as a vcf here:

http://ftp.ensembl.org/pub/current_variation/vcf/homo_sapiens/

CSHL-HAPMAP-HAPMAP-CHB.vcf.gz
CSHL-HAPMAP-HAPMAP-CHD.vcf.gz
CSHL-HAPMAP-HAPMAP-MEX.vcf.gz
CSHL-HAPMAP-HAPMAP-GIH.vcf.gz
CSHL-HAPMAP-HAPMAP-TSI.vcf.gz
CSHL-HAPMAP-HAPMAP-LWK.vcf.gz
CSHL-HAPMAP-HAPMAP-ASW.vcf.gz
CSHL-HAPMAP-HAPMAP-MKK.vcf.gz
CSHL-HAPMAP-HapMap-HCB.vcf.gz
CSHL-HAPMAP-HapMap-JPT.vcf.gz
CSHL-HAPMAP-HapMap-CEU.vcf.gz
CSHL-HAPMAP-HapMap-YRI.vcf.gz

Each SNP has an AF entry from which a multi-populations vcf with rsids, alleles and frequencies can be built, where the AFs are such as AF_CHB, AF_CHD, etc.

The 1000 genomes project populations documentation describes more than these 12 populations, but I haven't seen equivalent population vcfs with AFs built from the individuals within the population for the remainder of these, apart from the original HapMap 12 marked with a * below:

###Populations and codes

*         CHB   Han Chinese             Han Chinese in Beijing, China
*         JPT   Japanese                Japanese in Tokyo, Japan
         CHS    Southern Han Chinese    Han Chinese South
         CDX    Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV    Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
*         CHD   Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)

         CEU    CEPH                    Utah residents (CEPH) with Northern and Western European ancestry
*         TSI   Tuscan                  Toscani in Italia
         GBR    British                 British in England and Scotland
         FIN    Finnish                 Finnish in Finland
         IBS    Spanish                 Iberian populations in Spain

*         YRI   Yoruba                  Yoruba in Ibadan, Nigeria
*        LWK    Luhya                   Luhya in Webuye, Kenya
         GWD    Gambian                 Gambian in Western Division, The Gambia
         MSL    Mende                   Mende in Sierra Leone
         ESN    Esan                    Esan in Nigeria

*         ASW   African-American SW     African Ancestry in Southwest US
         ACB    African-Caribbean       African Caribbean in Barbados
         MXL    Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR    Puerto Rican            Puerto Rican in Puerto Rico
         CLM    Colombian               Colombian in Medellin, Colombia
         PEL    Peruvian                Peruvian in Lima, Peru

*         GIH   Gujarati                Gujarati Indian in Houston, TX
         PJL    Punjabi                 Punjabi in Lahore, Pakistan
         BEB    Bengali                 Bengali in Bangladesh
         STU    Sri Lankan              Sri Lankan Tamil in the UK
         ITU    Indian                  Indian Telugu in the UK

I presume there is more public data nowadays that can be accessed to build a more complete vcf files of allele frequencies per population, with as many populations as possible to use as a reference dataset with new data.

Thanks in advance.

The best set I'm aware of can be downloaded from https://gnomad.broadinstitute.org .


Random GWAS data generator I was wondering if there is a tool/script/program that randomly generates GWAS data. The purpose of such a tool would be to use it for educational purposes. So you generate some random .ped, .map and other filetypes. Various distributions like the MAF of all the SNPs, X-Chromosome classificaiton error rates and so on could be set before generating the data. General stuff one can then filter out with plink or similar software to increase the quality of the data.
Then a student would analyze these random files in plink and I already know what the outcome should be, therefore I can see right away if he or she did it right.
Does anyone know if such a program already exists? Or is there a somewhat easy way to program this myself?
Maybe it works if I use existing gene-data from NCBI where I use the given MAF to generate n variations the alleles based on that? Just some random SNPs from dbSNP?

Agreed with the above I would use an existing dataset from NCBIs dbSNP here
Random doesn't work well because HW signal will be meaningless as soon as the Jukes-Cantor distance is exceeded (which it will) and plink needs homologues to generate any signal either in 'reasonable' panmixia or clonality. The result will be panmictic beyond even the most panmictic organism ever known. I don't mind having a simulation of 100% clonalty, because 100% clonality does occur, either through inbreeding or binary division, but this simulation would generate the complete opposite result.


What is the significance of the effect size on WOMAC scores for TKA? https://www.medicaljournals.se/jrm/content/html/10.2340/16501977-2210

I understand the difference between statistical significance and effect size, but I am having trouble interpreting the data on Table III in the article linked above. I am evaluating the WOMAC scores. My first question is can effect size confidence intervals cross zero and still be valid? My second question is about the negative and positive values of the Hedges'g. The first study scores the WOMAC in a way where the higher the score, the better the outcome. The second two studies score the WOMAC in a way where the lower the score, the better the outcome. Does this mean that in the first study a positive Hedges'g indicate larger effect for the treatment group and a negative Hedges' g for the other two studies indicate larger effect for the treatment group?

Hedge's g measures the difference between two groups here experimental and  control in terms of combined (weighted) standard deviation. A higher g implies greater difference and lower g means lesser difference between  two groups.
It depends  on the study how the scoring of questions is proceeded with : the effect-size will be same. For example, answer to question item may be scored 1 to 10 depending on whether question item is written in positive  sense (the apple is good for health) or negative sense (the apple is not good for health). Positive Hedge's g does not indicate larger effect for the treatment group. THE negative Hedges' g for the other two studies indicates positive effect for the treatment group if the difference has been computed by  deducting experimental group mean from the control group mean (which is less than exp. Group Mean). Overall, I find that a knowledge of questionnaire scaling and scoring of scales and question items could be very helpful in general and in doing a meta-analysis effectively.


Add new tool to galaxy I am trying to include a new tool to the galaxy main menu, following these instructions.
tool_conf.xml doesn't exists only tool_conf.xml.main and tool_conf.xml.sample and none of them seem to be working. I make changes and nothing happens when I run galaxy.

What am I missing?

Thank you.

The .main and .sample file files usually come with default installations and may have been removed by whoever setup this Galaxy instance. You can look at the originals online and just base your own file on those: https://github.com/galaxyproject/galaxy/blob/dev/lib/galaxy/config/sample/tool_conf.xml.sample


Appropriate value to replace outliers in clinical studies I am working on clinical study where I am focusing on binary classification problem.

There are few variables which has outliers.

I am trying to know which would be the right value to replace outliers in a clinical study?

Is it mean, median min,max or mode of a column of interest?

Currently I use min and max for values below and above 1.5*IQR threshold

I use IQR approach to identify outliers but would like to know how it is done usually in clinical studies?

Can anyone help me with this?

In my experience, it is more common to use all clinical data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data.

If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.


Strategies to disentangle the effect of two genes on phenotype I am analyzing the effect of one putative tumor suppressor (and its loss) on a specific cancer phenotype. This putative tumor suppressor is known in literature to have an effect on the growth of tumors in mouse models.

However, in the same region of this gene there's a second well-known tumor suppressor - actually a very powerful one. By looking at the data from cBio portal, there's a very strong co-occurrence of loss on both tumor suppressors in a pan-cancer analysis, that is they are almost always lost together.

Therefore, I would like to disentangle the effect of the two genes, and to analyze only the effect of my putative candidate on the phenotype. The simplest approach is to use only the samples that have a single loss (only the loss of my putative candidate, but not the loss of the confounding tumor suppressor). However, only ~10% of all samples shows such behavior, and this 10% is distributed across all cancer types, decreasing the statistical power of the analysis - since I would like to test each cancer type separately to avoid mixing incompatible data.

Which kind of strategy/analysis can be used in this case?

Maybe you could fit a mixed model to predict phenotype, using the two genes as fixed effects and cancer type as a random effect, and ask the model for the separate effects of each gene, and the interaction between them.


How can I classify the 3 clades(S, G, V) of the coronavirus without using protein data? On GISAID they classified the coronavirus using 4 clades(S, G, V, Other).
I downloaded around 1,000 complete genomes of the coronavirus from GISAID and I would like to classify each one as belonging to one of the 4 clades(S, G, V, Other).
On the top left of the image below is the classification of the 3 clades(S, G, V).
More specifically, the three strains were classified based on three specific variants(This answer was taken from my previous question):

Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the  gene's protein product   to a valine (V).

My problem is that using this method of classification for the coronavirus genome I am seeking information on the protein products of the genes in the genome.
So to be more specific here are my concise questions:

How can I classify a genome as "Strain S" using only the complete genome assembly, without any information about its proteins?
How can I classify a genome as "Strain G" using only the complete genome assembly, without any information about its proteins?
How can I classify a genome as "Strain V" using only the complete genome assembly, without any information about its proteins?

This picture was taken from this PDF on GISAID

It seems to be explained right there in the image you posted:

So, the three strains were classified based on three specific variants:

Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the  gene's protein product   to a valine (V).