The algorithm behind makeblastdb? What is the algorithm/concept behind‚makeblastdb?

I know makeblastdb produces a .phr, .pin, and a .psq file, that represents an index or a data structure that is better optimised for searching a given query.

  What is exactly stored in these files?

Seems like .phr and .pin are storing some metadata (based on the size) and .psq the actual index.
It can't be the table of all k-mers since the command doesn't take a value of k? (or maybe it is fixed to 3)

I could not find any direct discussion of the algorithm since google is flooded with people asking how to build such a database...
Maybe someone can give me a link or hint?

From looking at this publication "Database indexing for production MegaBLAST searches" which although deals with MegaBLAST and it's indexing method, it would appear that this was rolled into the NCBI C++ toolkit as the command makembindex so I'm inclined to think the current implementation must have been derived from this strategy. Looking at the relevant part of the current NCBI C++ toolkit I find no reference to makembindex but only makeblastdb which I assume has superseded the older command: blastdb Directory Reference

Ulimtately the true details of the index generation lie in the C++ source code for makeblastdb can be found here.  However following on from that I think the further details of the index generation may be in one of the following files which are includes for the main .cpp file:

seqdb.hpp

build_db.hpp

seqdbcommon.hpp

writedb.hpp

It would appear from seqdb.hpp that:

The index files have extension .pin or .nin and sequence files .psq or .nsq,

However it's not at all clear how the current implementation differs to the 2007 reference or where the .phr file is generated from!  Some degree of reverse engineering the C++ source code will be needed!


Finding motifs: fasta file with 10,000 sequences I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that "...the object of type 'FastaIterator' has no len()...". The script however worked when only a single sequence was in the fasta file.
In essence, I need some help figuring out how to:

Prepare an instances object(is this correct?) with microsatellite
motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10
and (AG)10;
Parse my 10,000 long sequence file using this object(?); and
Print an output specifying which of my sequences contain microsatellite motifs.

Appreciate any help I can get on this.
#Seq Utils to find Motifs
>>> pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
>>> sequence = Seq(a)
>>> results = SeqUtils.nt_search(str(a),pattern)
>>> print(results)
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
 56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines

#BrutForce to find Motifs in a fasta file containing 10,000 sequences
>>> def BruteForce(s, t):
...     occurrences = []
...     for i in range(len(s)-len(t)+1):
...         match = True
...         for j in range(len(t)):
...             if s[i+j] != t[j]:
...                 match = False
...                 break
...         if match:
...             occurrences.append(i)
...
...     print(occurrences)
...
...

>>> t = "TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG"
>>> len(t)
34
>>> s = SeqIO.parse("group23.fas", "fasta")

>>> BruteForce(s,t)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "<console>", line 3, in BruteForce
TypeError: object of type 'FastaIterator' has no len()

#BrutForce to find Motifs in a fasta file containing a single sequence
>>> s = SeqIO.read("group23v.fas", "fasta")
>>> BruteForce(s,t)
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]

Bio.SeqIO.parse() returns a SeqRecord iterator. To access the sequence of this sequence record, you need to use the .seq attribute so you should update your s with s.seq.


Retrieve all genes for a specific gene ontology I'm having a bit of confusion about how can I retrieve all the genes that are associated to a specific gene ontology.

Lets say I need to find all the genes that are associated with this GO entry:

  GO:0030098 lymphocyte differentiation

This is the graph representing the parents of my GO. In the same page I find the Child Terms that represent the children of my term.

Now, I can download the goa file (related to the organism I'm interested into) containing the GO associations and grep my term to extract a list of associated genes but, as far as I can understand, this list does not represent the whole graph that actually describe my GO term.

Should I also get the genes in the childred terms of my term? Or in the parents terms? If yes, then again I have to go back in the graph and find the parents/children of the found parents/children and so on, till the whole graph is visited.

Is this the correct approach? Or should I only focus on those genes related to my GO:0030098, and that's it?

I found this:
Retrieving a list of human genes having GO associations
and this but I think they don't address my question since they do not describe the parent-children association.

Also I found this, on the GO FAQ, but what I need to do is quite the opposite.

Summarising: I need to find the genes associated to a specific GO term. Should I also look for its parents or children or not?

Any clues?

You can use BioMart to filter by your GO term and get the genes as attributes. BioMart is ontology-aware so will pull out all genes associated with your term and with any of its child terms.

There is no need to look up the GO child terms as BioMart already deals with this. A term may have more than one parent term, but if it's associated with a gene then everything about that term is associated with the gene. The term is not associated with the gene via parent X or parent Y, it means that both are relevant. So if the child term is associated with the gene, then the parent term is by proxy.

If you're not already familiar with BioMart, there's a help video.

Here's a BioMart query that gets you exactly what you need in human using the web interface

If you want to work with the R interface to biomaRt, here's some documentation and this is the identical query to above:

library(biomaRt)
ensembl=useEnsembl(biomart="ensembl", dataset='hsapiens_gene_ensembl')
getBM(attributes = c('ensembl_gene_id','external_gene_name', 'go_parent_term', 'name_1006'),
      filters = 'go',
      values = 'GO:0030098',
      mart = ensembl)


Power calculation for GWAS/EWAS I want to investigate, how much sample size i needed to obtain 80% power for GWAS/EWAS studies. Phynotype trait is discrete (not case/control) for human disease.
I wonder, does anyone has came across this? Really appreciate any suggestion or any tool/software which can answer this?

Interesting question I would recommend the effective population size $N_{e}$,  and error bars around this. Ne can be calculated via $F_{st}$ and is perhaps the most important parameter to determine the genetic population size.
Under the Fisher definition $N_{e}$ is the number of breeding individuals in the population.
The central issue with the approach you are suggesting is that population size does not equate to genetic population size. Thus the population could be huge, but if it is genetically homogeneous it will appear very small. In contrast the number of individuals could be very small, but the genetic variation between between them very large so $N_{e}$ would be huge.


Parsing specific features from Genbank by label? I'm trying to parse a Genbank file to find a specific feature. I can pull it out if I know the feature type (e.g. repeat_region). E.g. if I'm looking for this feature:

 repeat_region   5623..5756
                 /label=5' ITR
                 /note="5' ITR"

I know that I can find it using:

for feature in reference.features:
    if feature.type == "repeat_region":
        print(feature.location)

But I don't trust that it will always be a repeat_region. Instead, I'd like to look for it by label (5' ITR). I can seem to find a way to parse that from the feature object. Any suggestions?

The /label=5'ITR is called a qualifier. You can look for those qualifiers:

for feature in reference.features:
    for k, v in feature.qualifiers.items():
        if k == "label" and v == ["5' ITR"]:  # or "5' ITR" in v
            print(feature.location)


How to check what software has been used to generate bam file? I have several bam files from unknown origin. I'd like to know what software was used to obtain those bam files, e.g. bwa-mem. Is it possible to somehow check it?

If you're lucky the information is stored in the BAM header in the lines starting with @PG which stores the command lines run on the BAM, use samtools view -H your.bam to display the header.


Why is Mann-Whitney U-Test over an ExpressionSet yielding different results than original paper I am trying to replicate the results of this paper roughly guided by this pipeline.
Basically, we are trying to detect a differential expression between mesenchymal and epithelial cells under cisplatin treatment vs control.
However, I am facing a problem when trying to calculate the p-value for the Mann-Whitney U-test. This is the part of the paper I am referring to:

Robust Multichip Average (RMA) normalization was performed at the
transcript level on the results from the Affymetrix Human Gene ST1.0
arrays using Affymetrix Power Tool 1.14.4.1 for all 46 sham- or
cisplatin-treated ovarian cancer cell lines. The normalized data were
subsequently standardized using ComBat71 to remove the batch effect.
In this experiment, the cisplatin treatment assay was performed in
triplicate on 20 cell lines, while single assays (without replicate)
were performed on the remaining 26 cell lines. Taking advantage of the
triplicate data, potentially fragile probes with strong variations (an
s.d. of >0.2) within the triplicates were removed, decreasing the
probe number from 33 297 to 21 329. To perform a fair comparison, the
triplicate data were then log-averaged into one value so that one
result for each cell line could be used in the following analyses. To
detect differential responses to cisplatin between epithelial- and
mesenchymal-like cell lines, the transcriptomic responses to cisplatin
were computed by subtracting the gene expression value of control
(cisplatin untreated) cells from that of cisplatin-treated cells.
Mann–Whitney U-test (P<0.01 as a cut-off value) was subsequently used
to detect the differential transcriptomic responses between the
expression changes by cisplatin treatment in epithelial-like cell
lines with those in mesenchymal-like cell lines (Supplementary Table
3).

My p-values for a given gene are astronomically different to that of the paper...
Here is my code:
expr_f is the expression values matrix from the Expression set (log-transformed).
SDRF_f is the pData from the same expression set.
I am feeding into the wilcox.test function a vector of expression from epithelial cells (epi_values) and mesencymal cells (mes_values) for the same probe (8143663).
 genelist=rownames(expr_f)
epi_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Epithelial-like"))
mes_samples = rownames(subset(SDRF_f, SDRF_f[,"Classification"]=="Mesenchymal-like"))

epi_values = expr_f["8143663", epi_samples]
mes_values = expr_f["8143663", mes_samples]

wilcox.test(mes_values, epi_values, exact=TRUE)

In the paper, their p-value for this probe is 0.000372259258165627 while mine is 0.8785. The difference makes no sense at all to me.
I don't know if I am not quite getting the concept of the u-test, if I am using the function wrong, or if it's something else. Any input is appreciated

I think you may need to check your pipeline in this regard. I ran the analysis and replicated as closely as I could their M&M (I didn't filter the "fragile probes" since you anyway care only about a specific probe). In any case, I get something relatively close to their result (P=0.004059), which is still one order of magnitude off, but without more detailed M&M, this is probably as close as I can get to their analysis. Here's what I did (if you prefer to read the notebook html, that is here)
library(Biobase)
library(GEOquery)
library(sva) # ComBat
library(oligo)
library(googlesheets4)
library(stringr)
library(dplyr)

#' Get the GEO dataset mostly for the phenoData annotation
gset <- getGEO("GSE47856", GSEMatrix =TRUE, getGPL=FALSE)
adf <- Biobase::phenoData(gset[[1]])
sampleNames(adf)

#' Read RAW cel files from the GEO data supplementary rar archive
#' https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE47856&format=file
#' In my case extracted to ~/Downloads/cel
#' Ensure that we read them in the same order as the phenoData object
celfiles <- sapply(sampleNames(adf), function(f) {
  pat <- paste0(f, ".*.CEL.gz$")
  dir("~/Downloads/cel", pattern = pat, full.names = TRUE)
})
raw_data <- oligo::read.celfiles(celfiles,
                                 phenoData = adf,
                                 experimentData = experimentData(gset[[1]]))

#' Get the raw expression and RMA on transcript clusters
rma_expr <- oligo::rma(raw_data, target = "core")
rma_expr_assay <- assayData(rma_expr)$exprs

#' Use ComBat for batch effect removal
batch_data <- phenoData(rma_expr)[["batch:ch1"]]
cb_expr <- sva::ComBat(rma_expr_assay, batch_data)

#' Get supplementary table 1 from the paper (which I copied on this Google Sheet)
supp_data <- read_sheet("https://docs.google.com/spreadsheets/d/142__1O35zvlttejvz62XQX33uLRJhqIDvb604smVSRo")

#' Match the 'cell line' phenotype from the phenoData to the sample info of
#' the supplementary table and extract the cell classification. Also remove
#' all white space and convert everything to lower case for easier matching.
#' Some cell lines also have a/b attached, so I'm dropping that to match them
#' up. Mayb not be correct to do that here and treat all a/b the same (?)
cell_line <- str_replace(
  str_replace(
  str_remove(
    str_to_lower(adf[["cell line:ch1"]]),
    " +"),
  "^tay$", "taya"),
  "dov13[ab]", "dov13"
)

sample_info <- str_remove(str_to_lower(supp_data$Sample), " +")
classification <- supp_data$Classification[match(cell_line, sample_info)]

#' Get the treament from the phenoData object
treatment <- adf[["treatment:ch1"]]

#' We'll need the geometric mean later (R has no built-in)
gm_mean = function(x, na.rm=TRUE){
  exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}

#' Set up a tibble (==data frame) with all that info for your probe
tibble(Expression = cb_expr["8143663", ],
       Line = cell_line,
       Treatment = treatment,
       Classification = classification) %>%
  # Summarise replicated data with the geometric mean
  group_by(Line, Treatment, Classification) %>%
  summarise(Emean = gm_mean(Expression)) %>%
  # Make sure data is in order
  arrange(Line, Treatment, Classification)%>%
  # Calculate the difference between the cell types
  group_by(Line, Classification) %>%
  summarise(EmeanDiff = diff(Emean)) -> xdata

#' Test the difference
xdata <- split(xdata$EmeanDiff, xdata$Classification)
wilcox.test(x = xdata$`Epithelial-like`, y = xdata$`Mesenchymal-like`,
            exact = TRUE)

#' Session info
sessionInfo()
```


Finding new drug-protein targets/interactions using a known one as base Given a known interaction of a molecule with some proteases (having X-rays images too) is it possible to predict interactions with other proteins either using their 3D structure or not?

If anyone has already tried such techniques I would like to hear a little brief on that subject. Which are the possible techniques, computational power needed and so on.

If you have information on which part of the protein are interacting, and its sequence, you can search for homologous or similar sequences in other proteins. This would give you a list of potential targets and probably a similarity score which you can also use. Then you would have something to run simulations on.

And I think that with 3D structure you should be able to run simulations in silico in this case (eventhough it is probably more common to look for the ligands to the protein, not proteins to the ligand).

In silico simulations are computationally intensive, but I´ve also seen some being run on a laptop, so you don´t necessarily need access to HPC-clusters, or supercomputer.


removing mitochondrial read and unassembled "random" from multiple bam files This is the one i used for a single bam file to filter its mitochondrial as well as unassembled read

samtools view -h SRR2920466.bam | awk '{if($3!= "chrM" && $3 != "chrUn"){print $0}}' | samtools view -Shb - > this.filter.bam

Is there any way i can do it on multiple files with parallel as i do to index and sort bam files using samtools .

So i tried this ,which is of course not working as its not able to pipe properly because it shows me error regarding additional fonts

ls *.bam | parallel "samtools view -h {} | awk '{if($3 != "chrM" && $3 != "chrUn"){print $0}}'| samtools view -Shb {.}.bam"

Any help with my above codes? As i get it there is wrong with my pipe .

After updating the samtools to version 1.9 i ran the small script

i get this error

awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string

Your awk command is inside double quotes ("") which means that all of the awk field variables (e.g. $3) are being interpreted by the shell. Since $3 is not set in the shell, your awk is trying to compare an empty string and that's why you get the error:

awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error

Note how the error message shows if( != chrM and not if($3 != chrM. Since you need to use double quotes here because you need both an outer (for parallel) and an inner (for awk) quote, and single quotes can't be nested, you will have to escape the awk variables:

ls *.bam | parallel "samtools view -h {} |
    awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb {.}.bam"


FeaturePlot from Seurat: change its title FeaturePlot is a function in Seurat package.

And in the vignette it is written that if we specify parameter do.return = TRUE it should return ggplot2 object. It is not working. My goal here is just to change the title of the plot. In case of violin plot I can do the following:

VlnPlot(object = seurat_object, features.plot = id, do.return = TRUE) + labs(title = endothelial_symbols[1])

But with FeaturePlot similar code fails to work:

FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                          reduction.use = "tsne", do.return = TRUE) + labs(title = endothelial_symbols[1])

Giving the error:

  Error in FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey",  :
    non-numeric argument to binary operator

Any suggestions would be greatly appreciated.

I tried with some data that I have and this is working for me:

p <- FeaturePlot(object = seurat_object, features.plot = id, cols.use = c("grey", "blue"),
                      reduction.use = "tsne", do.return = TRUE)

lapply(p, function(x){x + labs(title = endothelial_symbols[1])})

I think it is because FeaturePlot returns several ggplot objects according to the ids you put. This is returned as a list of plots so you have to iterate each one for labeling.


The visualisation of a list of genes on URD object Seurat R package has some functions like FeaturePlot, DimPlot and DoHeatmap by which we can plot the expression of a list of genes on cell clusters. I am working with URD that likely does not have such options or I can not find that. URD has plotDot but would be clear only for a few number of genes (please look at the picture). Could somebody please help me in plotting the expression pattern of a list of genes in URD? As I am not good in R at all to write my own function or adopt Seurat functions to URD

For example Seurat DoHeatmap gives such a picture by which we can see these set of genes are being expressed on which group of cells

DoHeatmap(object ,genes.use=a vector of genes)

Or we can plot a set of genes on seurat feature lot to they are being expressed more on which cluster

FeaturePlot(object = pbmc_small, features.plot = 'PC1')

In URD only by

plotDot(object.6s.mnn, genes = a vector of genes, clustering="Infomap-60")

We can plot a set of genes but as you are being in the first picture that would not be clear much

My question is how I can plot a vector of genes on URD object

for example Seurat DimPlot function can visualise the expression pattern of a group of genes on cell clusters but when I am trying to do the same with URD, returns this error

> plotDim(axial.6somite, c("DDB_G0291149","DDB_G0291029"), plot.title="SHHA (Notochord marker)")
Error in names(data) <- cells.use :
  invalid to use names()<- on an S4 object of class 'dgCMatrix'
In addition: Warning messages:
1: In if (label.type == "meta" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
2: In if (label.type == "group" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
3: In if (label.type == "sig" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
4: In if (label.type == "gene" | (label.type == "search" & label %in%  :
  the condition has length > 1 and only the first element will be used
>
>

Actually only one gene can be visualised by this function on URD

I just added URD clustering labels for each cel to seurat object and then I did DoHeatmap by using group.by="URD clustering labels"


Identification of starting point for the synthesis of a protein identify where is the starting point for the synthesis of a protein in the DNA
For indentify a starting point in RNA or DNA my function search all ATG or M in the sequece and report the linear position of start, than the function search TAA and TAG for stop codon and report the end position.
For the end codon what I have to search? only TAA and TAG or TGA and TGG too?
How is possible some protein begin with "F" as in "FFRE..." and not "M..." ?
Than it result more starting point solutions respect to protein found in GeoBank, what are the add restriction for limit that number?
Are there easy to search sequences that make sure that is a start of the protein code? as seems to me "AGG AGG" but is not clear where they should be in the array and not seems general seen proteins.
Can the synthesis of protein be affected by temperature? What other factors come into play?
Passing my function to DNA of one virus, the first test, the function identifies
in the end of virus many little proteins of 2 3 4 5 10 protein letters. These little protein pieces are really sintetized? or there is a min lenght for sintetized proteins?

You seem to be interested in identifying potenital coding regions in a DNA sequence. These regions are called open reading frames (ORF) because they are open to translation.
As detailed in the linked Wikipedia article, an ORF usually starts with an ATG and ends with one of the stop codons (TAA, TGA or TAG in DNA space). However, there are many ORFs that are not translated at all, particularly short ones.
Since finding ORFs is an old and established task in bioinformatics there are many algorithms availabe to do this. I would suggest using such an existing tool because they often have additional features and bugfixes. There is a list in the link above.
Estimating the stability and likelihood of a particular amino acid sequence is much more complicated. There are many tools for protein structure prediction that may be worth having a look at.


Does SBOL require representing intermediate products like mRNA? Do I need to include intermediate gene products like mRNA in order to express a gene regulatory relationship in SBOL?

For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR. Do I need to represent the TetR transcript mRNA as well as the TetR protein?

As far as I understand, the purpose of SBOL is not to create a comprehensive model, but to communicate a design.  So the question you have to ask is "is the mRNA part of my design, or is it an implicit part of the biology which does not help to convey what I'm trying to build?".

I think your summary is a giveaway:

  For example, let's say that I'm representing transcriptional regulation of the pTet promoter by means of a coding sequence that expresses TetR.

Notice that the mRNA was not mentioned?  I also daresay that if you drew this design as a diagram, the mRNA would be omitted - another indication that it is not helpful in communicating what you designed.


Make a fasta out of mapped reads without taking into account the reference I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide.
Is there a way to do this ? Thanks

vSNP outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  vsnp_merge_vcf_into_fasta.py.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".

vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

Use *zc.vcf output to merge into FASTA

vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf


As you surmised, kallisto requires a transcriptome fasta file, not a genome fasta file.


Writing a python module to run Infernal on a list of fasta files I am very new to using python3 and am hoping to use it to reduce the time I need to spend running Infernal to look for a variety of RNAs in a large collection of bacterial genomes. I am hoping to write a short module in python that will let me sequentially run infernal using a single covariance model on a list of genome files.

At the moment I am trying to use subprocess.run() or subprocess.Popen() to do so.

import subprocess
import os
import glob

genomes = list(glob.glob(os.path.join("/path/to/genomes/", "*.fa")))

for x in files:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", "genomes"], shell=False)

my issue is I do not know a way to have "genomes" passed as the components of the list within the command.

Is this a way off approach to trying to automate this? If so is there a better way to do it or is there a way to simply pass each file in the list to the cmsearch command?

I think what you want to do is:

import subprocess
import os
import glob

genomes = glob.glob(os.path.join("/path/to/genomes/", "*.fa"))

for gfile in genomes:
    subprocess.Popen(args=["cmsearch", "/path/to/covariancemodel/", gfile], shell=False)

each pass through the for loop will take a one of the genome files and pass it into the cmsearch command.

Note:

1) glob.glob returns a list so no need to cast it again - https://docs.python.org/2/library/glob.html


Since every human has a different DNA (different combinations of C, G, A, T) what does it mean to have the genome done? I'm confused about the difference between genome and DNA. Is it correct to say that the same type of bacteria has the same DNA? But my understanding is that it is not correct to say that the same type of human has the same DNA, since every human has a different DNA. What am I missing here?

All humans have some differences in their DNA, but there's far more that is shared. On average the difference between humans is only about one thousandth of their full DNA, which means we're about 99.9% the same. These differences aren't distributed fully randomly, but are often because of specific gene alternatives. (Random mutations do occur, but they are also often fatal, so the random mutations we see of living adults are much more restrictive than all the random mutations that occur.)

To identify the human genome is to study many people's DNA and to label the parts that are shared between everyone, the parts with two or three variations, and the parts with even more variations. Even though every person has different DNA, we can still say they fit the pattern, just like every T-shirt may be unique but they all fit the T-shirt pattern and not the trousers pattern.


I have a list of protein targets. How to get structural data (atomic coordinate) from PDB? I have a list of around 500 protein targets (human). I only have their protein symbols. As an example:

EGFR
AR
FOXA1
CXCR3
...

I want to get structural data (atomic coordinate) from PDB (Protein Data Bank) at once, either programmatically (Python or R) or with a web api. Do you know a clear step by step tutorial for this task or code example?

If that is not available, at least I need to find the PDB IDs of my list of protein symbols at once.

thanks

There is not a simple way and there are factor that you also need to consider before starting.

Are you using gene names or gene accession numbers/id (NCBI/Uniprot)?
If not, are the names the current standard ones (e.g. scraped names from old papers is a definite no)
PDB residue index is rarely the same as protein index.
Do you want a perfect match or do you want say 80% homology?
Which PDBs do you want for a protein?

Personal example

For a website I have developed I have a route that gets all the PDBs for a given gene. (enter link description here). The way I did it is that I parsed Uniprot using a custom parser and made a dictionary/DB for each species with many synonyms for a gene going to a Uniprot ID and then from there I have another dictionary/DB that goes from Uniprot to PDBs if any. I went for Uniprot as I needed other data too (e.g. features) to present to the user. Including a renumbering option.
All this is way overkill for what you need. However, if you specify the species I can easily share a JSON file online.

SIFTs

If point 3 is important and you have Uniprot IDs I would suggest going with SIFTs. It has the offsets and has mapping for Uniprot IDs (pdb_chain_uniprot.tsv).

Gene names

In terms of mapping gene names to NCBI/Uniprot/Ensembl IDs, there are mapping tables everywhere. But if you may have obscure synonyms, using the NCBI API is a good option.
Uniprot data is better for protein and has PDB codes and chains, but the API gives a XML reply which is a bit annoying (I wrote a parser for it, but I think there was going to be one in Biopython).

Homologues

If you would like homologues, say your organism is human and you are okay with mouse, you can either use the threaded SwissModel models or use the blast query in RCSB PDB (not the PDBe, which has better data normally) or NCBI set to PDB database.
Swissmodel actually contains also PDB structures, so that may be a good starting point if your organism is one of those where precompiled exists.

Parts

Some proteins, especially mammalian are formed by separate domains linked together, especially the docking scaffolds. In which case you may care for a specific model within the gene, hence the indexing issue and the potential getting features.
If you want a single PDB structure for each gene and do not care about the technical I have some Pymol-using py3 code to fuse models together with the N and C termini aligned on line. But I strongly discourage that and would argue for protein feature based routes.


Separate multiple sequence into different file, one sequence per file I have a file contain multiple sequence, and I want to separate them by "gene:" into different file.
example:

example.fa

>KQK21959
ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

>KQK14546
ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

how to get results like this? one sequence per file, and after ">" is gene id, not cds id

KQK21959.fa

    >KQK21959
    ATGGTGCTCGACCCCGTCGTCGGGGGCTCGGCCTCCCGCAAGGTCGAGCCCCTAACCCTA

KQK14546.fa

    >KQK14546
    ATGCCTCCGCCGGACGCGGCCACCTACGCTCCCACGCTGCCAATGGCGACGTACCTCCCA

There's an app for that! The exonerate package (see here) includes a nice fastaexplode command that does exactly what you're looking for and will be much, much faster than anything we give you here since it is written in C and not a scripting language. Just install exonerate (you can do that with sudo apt install exonerate on Debian-based systems) and then:

fastaexplode file.fa

This will create these two files (file.fa is the original input file from your question):

$ ls
file.fa  KQK14546.fa  KQK21959.fa


Per Base Sequence Content in fastqc I have a question regarding "Per Base Sequence Content" plot for "fastqc":

In the fastqc documentation, it is written: "In a random library you would expect that there would be little to no difference between the different bases of a sequence run, so the lines in this plot should run parallel with each other."

But I don't understand why different bases in a read should follow the same pattern of allele frequency ("A/T/C/G"). I mean they are different positions in the genome and it is normal that each position has different allele.

I would appreciate it if someone could help me, please.

The positions you see in the x-axis of this plot are positions on your reads, not on the reference genome. Since each read comes from a random position on the genome, the frequency of A/T/C/G on each base position should reflect the base composition of the entire genome. Therefore, the percentage of each base should be approximately constant across the read.

Hope that helps!


Downloading SRA Files from AWS I want to download the original BAM files that the authors had uploaded to SRA. Normally, I would just use sam-dump, but the files are having issues that seem related to this issue. Since according to the entry, AWS S3 also hosts the original BAM files, I thought I could download these directly.

Example SRA Entry

NCBI documentation implies that I can't download this directly, but I can freely copy to other AWS locations within the region. To this end, I created my own S3 bucket (mm-mneuron) and am now trying to copy from the SRA bucket to mine.  Here's what I try:

import boto3
import botocore

s3 = boto3.resource('s3')

bam_file = {
  'Bucket': 'sra-pub-src-6',
  'Key': 'SRR5253957/RPI25_0.bam'
}

my_bucket = s3.Bucket('mm-mneuron')

my_bucket.copy(bam_file, 'RPI25_0.bam')

This fails with:

  botocore.exceptions.ClientError: An error occurred (403) when calling
  the HeadObject operation: Forbidden

That is, it sounds like I can't access the SRA bucket. I've tested downloading and uploading to my bucket, so I know I have write permissions. Not sure what else to try here.

How can I access the SRA data on S3?

A member of the SRA submission staff pointed out that using

prefetch --type all SRR5253957

will download the original files. In this case, it means running the above within an EC2 instance colocated with the S3 bucket (so, us-east-1) and having installed and configured SRA Toolkit to work from AWS (as per this documentation).

Unfortunately, the particular files I am concerned with are not currently accessible, but this should generically work in most situations.


Using `install.packages` with conda-managed R I have R installed and managed by conda (miniconda) on my MacBook Pro. The version of R I use most frequently (3.5.1) is installed on the base environment and I have other version-specific environments as well. I did this as it was easier than managing R on homebrew.

The downside I now face is that I cannot use install.packages(), BiocManager::install() or biocLite() to install most packages. devtools::install_github() seems to work fine for the few packages I have installed from GitHub.

I have to use conda install -c [conda-forge|bioconda] [r|bioconductor]-package_name to install my packages and it breaks my workflow and I need to step outside of R to install packages.

Am I doing anything wrong, or is this the only way packages can be managed when R in installed through conda?

A few sample packages that I needed to install using conda:

tidyverse
tximport
deseq2
apeglm
maftools
sva

Installing packages outside of an R session is arguably the point of using conda.  Package management isn't something that should be done within scripts, only the loading of packages. Doing so ensures that your environment setup is somewhat robust, which should be set before you perform actual calculations.

Using conda install ... is the best way to ensure environment stability here, and, thankfully, lots of people at Bioconductor have worked very hard to make sure bioconda packages fit with the conda framework.

If you need to find packages, you can search online at https://anaconda.org, or from the command line via conda search <pkg>.


How can I obtain FTP links to studies in ENA? How can I programmatically obtain ftp links to RNA seq fastq files in ENA? Here's an example of a link that I would be interested in obtaining:
ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR824/000/SRR8240860/SRR8240860_1.fastq.gz
In particular, is there some tool that, given the BioProject ID (here, PRJNA506829), would be able give me all ftp links for the runs in the project, or would I need to write a web scraper to do it?

Here are some command lines I used for that purpose in bash.
Simply prepare a text file containing each accession number (SRR/ERR) you want and create a for loop. Here I used prozilla to speed up downloads but you may use wget either.
for index in $(cat list_of_accessions) ; do

if [ ${#index} -eq 9 ]; then
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_1.fastq.gz
            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index}/${index}_2.fastq.gz
            else
            if [ ${#index} -eq 10 ]; then
                    proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_1.fastq.gz
                            proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/00${index:9:9}/${index}/${index}_2.fastq.gz
                    else
                    if [ ${#index} -eq 11 ]; then
                     proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_1.fastq.gz
                                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/0${index:9:10}/${index}/${index}_2.fastq.gz
                    else
                        proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_1.fastq.gz
                                                proz -k=6 --no-curses ftp.sra.ebi.ac.uk/vol1/fastq/${index:0:6}/${index:9:11}/${index}/${index}_2.fastq.gz
                        fi
                fi

            fi
done
```


Alignment to verify plasmid from forward and reverse Sanger sequencing reads I am analyzing Sanger sequencing data to verify whether or not the correct protein coding sequence was inserted during cloning.
For each sample, I have both forward and reverse Sanger sequencing reads that cover the coding sequence.
Would it be better if I first align the forward and reverse reads to get a consensus, and then align the consensus to the reference?
Or should I just align each read directly to the reference?

Align all reads from a sample at the same time IMO. This will give you all of the information at once, and it should not be difficult. It should be very obvious whenever there's a disagreement between reads or between reads and reference.
I would also recommend inspecting chromatograms directly at all sites of interest (such disagreements, insertion sites, possible errors or mutations).


Fastq: how can I check if they are from DNA or RNAseq data? I have (gave me) Illumina fastq files, which I want to use for variant calling, and I do not know if they are DNAseq or RNAseq data. How can I check this? I do not have any report or who to ask.
Many thanks

There are a few quick'n'dirty ways depending on the type of data.
In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC
which calcualtes the fraction of reads aligning to exon, intron, intergenic etc.
RNA-seq, if you use a standard aligner such as bowtie2 or bwa, and you have a higher species that splices its RNA, then you should see most reads aligning to exons, and quite a fraction of reads that span a splice junction being unmapped (because the exon-intron-exon gap is large and non splice-aware tools do not "bridge" this without tremendous drop in mapping quality so the reads goes unmapped).
DNA-seq (if whole genome) should have a somewhat even coverage across the genome on the global scale. Targeted assays such as ATAC-seq should give distinct peaks that you can easily see in a genome browser just by eye.
If it is exome-seq (so targeted DNA-seq via exon pulldown) then it might get tricky because also most reads (like in RNA-seq) should align to exons, but you should not have spliced reads so the 5' part of the read starting in exon1 and 3' part ending in exon2. Instead many reads should also cover the exon-intron boundary.
As said, I would simply align it and check on a browser, I could imagine one can make a sophisticated guess by just looking at the alignments by eye.
A different topic is reliability, because I personally would not touch a dataset with unknown origin, because that also means you have no idea about the protocols that were used for library prep and what the sample actually is, so I wonder which question your analysis is going to answer.


Can an alternate allele be more common than a reference allele? I'm trying to analyze data from a genotype-phenotype association study. The genotype data is stored in csv files that record for each variant site and each subject, the alleles that the subject has at that site. A quick online search reveals that 0 means the reference allele, 1 means the most common alternate allele, 2 means the 2nd most common alternate allele, etc.
But I noticed that at some variant sites, most (>90%) subjects have 1/1, i.e. two copies of the alternate allele. I asked around and some suggested I just swap the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5. My questions are:

Can there be a good reason to classify a less common allele as the reference allele? I thought that the reference allele was defined as the most common one in the population.
Is swapping the roles of reference and alternate alleles whenever I find that the alternate allele frequency > 0.5 a standard practice?

Thanks

Yes, depending on how the reference allele was determined. If, for example, you want to refer to hg38 as the reference genome, then it's important to use the hg38 allele as reference.

Without stating the reference allele convention? Not that I'm aware of. However many people refer to minor/major alleles rather than non-reference alleles, in which case it is assumed that the major allele is the most common allele.

As long as you state what your reference genome and/or population is, then it's fine to choose your own convention about reference alleles. But please be consistent and explicit about what you've done. One of my biggest struggles in GWAS has been working out the correct strand orientation for stated variants (important when merging datasets), and it doesn't help for complementary variants (i.e. A/T, C/G) when the ref/alt alleles are arbitrarily chosen.


Link to UCSC gtexGene track description page for selected gene symbol I have this link that comes up when I am in the browser window and click on the gene name, which leads me to description page about the gene (in this example it is BOK gene).

How can I generate that link? Below are some of my failed attempts. There is a UCSC FAQ page about the links, didn't help much, maybe I am reading it wrong.

Link that works, generated when clicked on gene on the website:

http://genome-euro.ucsc.edu/cgi-bin/hgc?hgsid=226968178_C6XmUjsjtt3Mm9fOmePo6VqRN6TV&c=chr2&l=241667282&r=242920780&o=242498135&t=242513546&g=gtexGene&i=BOK

Which gives me this page:

My failed attempts to generated custom link:

http://genome-euro.ucsc.edu/cgi-bin/hgTracks?db=hg19&position=chr2%3A241667283-242920780&gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&hgg_gene=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&i=BOK

http://genome.ucsc.edu/cgi-bin/gtexGene?org=human&db=hg19&g=gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/org=human&db=hg19&g=gtexGene&i=BOK

http://genome.ucsc.edu/cgi-bin/hgTracks?org=human&db=hg19&g=gtexGene&i=BOK

You never want to include the hgsid parameter in your url because that is the link to your current browsing session and all its specific contents. So if you send that link to someone else and then they start browsing around that link may not work correctly again in the future if you or the person you sent the link to modify your browsing session.

As to how to create links to hgc pages, the usually follow a somewhat general format:
http://genome-euro.ucsc.edu/cgi-bin/hgc?db=hg19&c=chr2&o=242498135&t=242513546&g=gtexGene&i=BOK

You need the db and g parameters to specify the assembly build and track, and then you need to provide the position in the form of the c,o and t parameters. The last thing you want is the i paramter to specify the item name. The reason the i parameter alone is not enough is because we do not require uniqueness on item names (there may be many copies of an item across the genome for instance).

You can obtain the c, o, t and i parameters by selecting the chrom, chromStart, chromEnd and name fields from the gtexGene table:

$ mysql --user=genome --host=genome-mysql.soe.ucsc.edu -A -P 3306 -e "select chrom, chromStart, chromEnd, name from gtexGene limit 1" hg38
+-------+------------+----------+---------+
| chrom | chromStart | chromEnd | name    |
+-------+------------+----------+---------+
| chr1  |      11868 |    13052 | DDX11L1 |
+-------+------------+----------+---------+


Create galaxy local toolshed using Planemo I tried to create galaxy local toolshed using planemo. I used command as follows,

create virtual env

virtualenv .venv; .venv/bin/activate

Install planemo

pip install planemo

To configuring a shed account

planemo config_init

Create a repository

planemo shed_init --name=seqtk_seq \
                  --owner=planemo \
                  --description=seqtk_seq \
                  --long_description="Tool that converts FASTQ
                                      to FASTA files using seqtk" \
                  --category="Fastq Manipulation"
planemo shed_create --shed_target local

This produces the following error:

Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.
(.venv) varuna@varuna-Inspiron-15-3567:~/git/galaxy_stable/planemo

--shed_target local will only work if you're running a Galaxy toolshed on your local computer, since it's reading the .shed.yaml file you made with planemo shed_init and trying to create a repository with that name or a local (in this case) server. Planemo itself won't start a Galaxy toolshed, since it's fairly unusual for people to want to do testing involving creating local repositories. What you probably want to do instead is use planemo test instead. Once the planemo testing passes then use the various shed* commands to upload your wrapper to the toolshed or testtoolshed.


What is outliers and how to determine outliers in the following phylogenetic tree? I an confused with the term 'outlier' and 'outgroup' in phylogeny. Could someone please explain with the example phylogenetic tree above which is the correct term? Thank you very much !!!

Distance- and ML-based algorithms using reversible models can't find the root of trees. A classical method to root a tree is to use an outgroup (not outlier), which is a species/sequence known to directly descend from the root of the tree. In your case, it is relatively easy: add a chicken or fish ortholog to your dataset and put the root on the chicken/fish branch.

Outgroup works well for single-copy genes used in species tree construction (e.g. ADH). However, it doesn't always work. The culprit is the "known" part. Say a gene has two copies A and B in vertebrates. Copy A was lost in rodent and copy B lost in primate. If you choose a chicken A gene as the outgroup, the correct tree should be ((primary-A,chicken-A),rodent-B). Without knowing the true history, you may forcefully put the root at the chicken-A branch and build a wrong tree ((primate-A,rodent-B),chicken-A).

There are a few other tree rooting methods. An easy approach is to put the root at the longest branch in the tree, assuming the presence of molecular clock. When the species tree is known, you can root a gene tree by minimizing the number gene duplication/loss events in the history. I generally prefer the latter approach when the relevant information is available.


how can i search in p-blast against UNIPROT(not only swissprot) i need to search for one protein against UNIPROT in BLAST for a selected organism.
How can I do it? (SWSISSPROT+TREMBL)
THANKS!!

I don't know of any online tool enabling a search against TrEMBL. You probably will have to set up your own local BLAST environment, roughly following these steps:

download blast+ - but there might be an easier way to install the software, depending on your operating system (e.g. on Ubuntu you can just sudo apt install ncbi-blast+)
download target databases
build local database

I hope this helps.


What are the algorithms for biometric identification in large database? I have been thinking about this problem and reading some related papers for sometime.

The problem can be described as follows:

Suppose that each biometric template is a $d-$dimensional vector. A database is composed of $n$ such templates ( $n$ could be large, like 1 million). Given a query template $q$, the task is to identify the "closest" template according to some distance metric.

Now, since the dimension $d$ and the database size $n$ could be both large, it is not quite efficient to linearly scan the whole database.

So my question is, are there any algorithms to efficiently fulfill the above task? What are the state of the art algorithms for biometric identification?

Through some literature review, I can think of two approaches. But both of them seem to have some limitations.

Building space partitioning index like R-tree. The problem is that such index does not suit for high dimensional data. The final search efficiency would be close to linear search.
Locality sensitive hash. The problem is that this method can only provide approximate searching result. For biometic identification, it is not guaranteed that the matched template could be found.

I think that a very nice example of what you are talking about saying the word bigger, is probably a microarray transciptional analysis of RNA extracted from a biological sample. Simplifying, in these analyis, each observation is analised on, up to 10000 genes expression levels of different mRNA coding genes. Thus, from these type of molecular biology experiments, we obtain matrices with few observations (d, the samples) respect to the dimensions observed (k, the target gene expression levels). Well, d << k. In this situation, in order to perform a classification-like task, and to reduce the problem of the overfitting related with dimensionality maledition, the ways that we might follow are several:

1) PCA: it is a good tool to visualize multidimensional datasets reducing dimensionality and noise level in the consolidated and standardized data. Problem: it is an unsupervised machine learning algorithm. Alone, it haven't the capacity of class discrimination.

2) PCA-LDA is a good combination, where PCA eliminates  hyperdimensionality and background in our data sheets, and than LDA (linear discriminant analysis) finds the directions onto the new subspace of characteristics, that are able to maximize variances between classes and decrease the variances within these. For this reason LDA is known as a kind of supervised machine learning technique. Furthermore, PCA and LDA are usually used together because LDA applied on high dimensional datasets could lead us toward overfitting problems. After the training of PCA-LDA pipeline, being LDA a classificator object trained and shaped on the principal component subspace, you can use this to examine new unknown samples or observation and predict its class by similarity. e.g. recently, Martin et al., on a Nature article talk about PCA-LDA approach to identify and than to predict cancer cells from those are normal, starting from ranman spectra of tissues.

3) PCA-LDA on transpose matrix: to use as better as possible the huge dimensionality, some authors used to transpose the initial matrix in order to force positioning of genes as rows and samples as columns. Then, the subsequent LDA is performed not on PC scores, but on loading factors of the eigenfactor matrix obtaining thus the list of the most important discriminant gene expression levels to segregate normal samples to cancer ones.


Software for taxonomic assignment? I have a couple of hundred bacterial sequences of 2-30 genes of interest each, recovered from metagenomics. None of them encode rRNA. Normally I'd just BLAST the one gene I already know to be reliable phylogenetic marker, where present, but the scale of my current dataset makes it a bit daunting.
Is there any software for taxonomic assignment that you could recommend? It's fine if it requires command line use. There are so many options listed online that I didn't even know where to start... Please and thank you.

You could take a look at kraken2 It should fit your purpose.


How to get Nanopore MinION fast5 from SRA I found some Nanopore MinION data on SRA, which I would like to investigate. I use sratoolkit for Illumina data all the time, but I am not sure how to get the fast5 file from the sra file.

What I have tried is to prefecth it, but then I don't know which command to use for .sra to .fast5 conversion. Not with fastq-dump, I assume?

prefetch -v SRR5286963 --max-size 50000000

When I was trying to upload MinION data to SRA, they wouldn't allow me to upload only the raw signal data (i.e. FAST5 files). Every file needed to contain a FASTQ sequence, and this sequence was extracted from the file and stored in the archive. Given this, I don't expect it will be possible to retrieve signal or event information from SRA. I also don't expect this will be possible in the future either, given that the SRA database expects every record to have a FASTQ-like sequence attached.

However, EBI/ENA do allow the upload of raw FAST5 files to the database (e.g. here), but there are no FAST5 reads for the specific single-cell B1a study you have indicated. If you're interested specifically in that study, then I recommend that you ask Miten Jain directly for access to those reads; he should have another method to give you access to them.

ENA is "ready" for ONT reads, but researchers need to upload their reads to ENA in order for that database to be useful.


Is it possible to identify cells that are expressing two or more genes in Seurat? I'm interested in looking into cells that are positive for two (or in some cases more) genes. I know I have some double positive just by looking at the FeaturePlot of those genes, but now I'm trying to figure it out how many are double positive and further percentage wise. Example:

50 cells are positive for Myh11
100 cells are positive for Acta2
30 cells are double positive
60% of Myh11+ cells are Acta2+
30% of Acta2+ cells are Myh11+

Is the data used to generated the FeaturePlot the scale.data? If so, I could just make a new dataset out of that, but I'm not sure if that is the case.

The counts stored in the Seurat object are: raw counts (seuratobject@raw.data), the log + normalized counts (seuratobject@data), and the scaled counts (seuratobject@scale.data). FeaturePlot() plots the log + normalized counts.

In order to identify double-positive cells, you need to identify cells that express a gene (i.e. positive for a gene) and that is not trivial: you need a cutoff value for each gene, and also to distinguish between true and dropout zero counts.

Provided you solve these, you can retrieve your cells by:

MYH11.cutoff <- 1
ACTA2.cutoff <- 1
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff))
length(which(FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))
length(which(FetchData(seuratobject, vars.all='MYH11') > MYH11.cutoff & FetchData(seuratobject, vars.all='ACTA2') > ACTA2.cutoff))


How to get the intersection of peaks after peak calling using MACS2? I have following 5 narrow peak files after peak calling.

K14_peaks.narrowPeak
K15_peaks.narrowPeak
K16_peaks.narrowPeak
K3_peaks.narrowPeak
K8_peaks.narrowPeak

I need to get the intersection of peaks into a new file (The peaks that represent all these files). I used following code but I got higher number of peaks not a lower number.

bedtools intersect -a K14_peaks.narrowPeak -b K15_peaks.narrowPeak K16_peaks.narrowPeak K3_peaks.narrowPeak K8_peaks.narrowPeak -wo > intersected.bed

Can someone help me?
Thank you.

bedtools intersect with multiple files after -b performs pairwise intersections between -a and each file in -b.
It does not intersect all files simultaneously.
See here, from the Bedtools documentation. The green section only contains intervals from -a that also intersect at least one interval in any -b file.

What I think you're looking for is bedtools multiinter (documentation here is unfortunately blank).
bedtools multiinter -i input [input...] will produce a BED file that covers every position listed in your input files, followed by how many and which of your input samples overlap that given position.

You can use this BED file and create something like an UpSet plot to visualize how many bases/peaks are shared between any combination of your samples.


Tag density plot form chip or atac seq data I m quoting this  paper , figure B where they write

"B) Violin plots for the normalized ATAC-seq tag density for all peaks within the indicated OC/CO groups at the indicated time points. Data were converted to a Z score based on the row-wise SD for each peak. ∗p < 0.01 Mann-Whitney U test."

My question is where is this tag density information coming from.

So far for chipseq data what I have done is

Aligned the data
Then all required data processing post alignment
Generated bed and bigwig files of respective samples required for comparison
Made tag-directory using homer.
Then peak-calling using homer,annotated peak using homer.

Now so far I am not able to figure out where is this tag-density information is .

Any help or suggestion would be really appreciated.

You can use Homer to retrieve tag densities, quoting from their manual:

  Calculating ChIP-Seq Tag Densities across different experiments

  annotatePeaks.pl is [a] useful program for cross-referencing data from
  multiple experiments.  In order to count the number of tags from
  different sequencing experiments, you must first create tag
  directories for each of these experiments.  Once created, tag counts
  from these directories in the vicinity of your peaks can be added by
  specifying "-d   ...".  You can
  specify as many tag directories as you like.  Tag totals for each
  directory will be placed in new columns starting on column 18.  For
  example:

  annotatePeaks.pl pu1peaks.txt mm8 -size 400 -d Macrophage-PU.1/
  Bcell-PU.1/ > output.txt


Which values from Seurat::FetchData function are to be used for correlation analysis between genes? I want to perform a correlation test between genes in on my single cell RNA seq data set. I perfomed the differential expression analysis using the Seurat version 2 package, after performing stages of normalisation, scaling, PCA, TSNE analyses and clustering. My correlation plot is far from what I expect and I want to know if the expression values obtained with Seurats V2's FetchData function (imputed, scaled, or raw) may have an effect on correlation analysis' outcome. This is what I did:

# load Seurat v2
library(Seurat)

# load corrplot package
library(corrplot)

# obtain expression values of the genes I want to test rom the Seurat object
expression_matrix <- FetchData(object, vars.all = c("gene1", "gene2", "gene3"))

expression_matrix <- round(expression_matrix , 1)

cor_expression_matrix <- cor(expression_matrix)

corrplot(cor_expression_matrix , is.corr = T, method = "square")

Which gene expression values should be best used for the analysis of correlation between genes: imputed, scaled, or raw?

Thank you in advance for your kind response.

You should use normalized data. The function has a parameter slot="data" that let's you select to access the normalized data.
Any reason for rounding to 1 decimal place here? I imagine that would have a substantial impact on the correlation you compute.


bowtie2 options when mapping stranded single end reads I'm trying to map RNA-Seq reads generated using the NEB Ultra Directional kit is first strand reversed protocal. I'm using bowtie-2.3.5_1 to align against the mus musculuous (mm10) genome.  The reads are single ended.

It seems that --fr and --rf are relevant only for paired end options. I tried these options as well as not specifying either (default).

E.g.

$ bowtie2 --fr \
        --threads 20 \
        --very-sensitive \
        -x /reference/mus_musculus/mm10/ucsc_assmebly/illumina_download/Sequence/Bowtie2Index/genome \
        -U trimmed_out/my_trimmed.fq.gz \
        -S bowtie_fr_1.sam

$ samtools sort -n -o bowtie_fr_sort_1.sam bowtie_def_1.sam

I then compared the output of bowtie_rf_sort.sam bowtie_default_sort.sam and bowtie_fr_sort.sam.  They are basically identical excluding the header.

Question:

Does strandedness matter in this alignment? Perhaps I should care about the --nofw/--norc options? Am I missing something here?

I repeat what user172818 says...don't waste your time on bowtie or tophat.  Use STAR or HISAT2, or a pseudoaligner like kallisto with a transcriptome.


Quast duplication ratio and mismatches percent When analyzing Quast results it seems that it doesn't calculate mismatches and indels in a useful way if the "Duplication ratio" is over 1.
For example, that's what I get for an assembly with two contigs covering almost the full genome

And that's what I get if I remove one of the contigs:

See how the mismatches and indels per 100kbp are almost double in the first case but that does not correspond to the reality.
Any better way to get these metrics? Would it be safe to divide error rates by the duplication ratio?

I took a quick look at the quast code that I believe is responsible here:
def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
    indels_info = IndelsInfo()
    genome_mapping = {}
    for chr_name, chr_len in reference_chromosomes.items():
        genome_mapping[chr_name] = [0] * (chr_len + 1)
    with open(used_snps_fpath, 'w') as used_snps_f:
        for chr_name, aligns in ref_aligns.items():
            for align in aligns:
                ref_pos, ctg_pos = align.s1, align.s2
                strand_direction = 1 if align.s2 < align.e2 else -1
                for op in parse_cs_tag(align.cigar):
                    if op.startswith(':'):
                        n_bases = int(op[1:])
                    else:
                        n_bases = len(op) - 1
                    if op.startswith('*'):
                        ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                        if ctg_nucl != 'N' and ref_nucl != 'N':
                            # MISMATCHES INCREMENTED HERE
                            indels_info.mismatches += 1
                            if qconfig.show_snps:
                                used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += 1
                        ctg_pos += 1 * strand_direction
                    elif op.startswith('+'):
                        # INDELS INCREMENTED HERE
                        indels_info.indels_list.append(n_bases)
                        indels_info.insertions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = '.', op[1:].upper()
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ctg_pos += n_bases * strand_direction
                    elif op.startswith('-'):
                        indels_info.indels_list.append(n_bases)
                        indels_info.deletions += n_bases
                        if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                            ref_nucl, ctg_nucl = op[1:].upper(), '.'
                            used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                        ref_pos += n_bases
                    else:
                        ref_pos += n_bases
                        ctg_pos += n_bases * strand_direction
                if align.s1 < align.e1:
                    for pos in range(align.s1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
                else:
                    for pos in range(align.s1, len(genome_mapping[align.ref])):
                        genome_mapping[align.ref][pos] = 1
                    for pos in range(1, align.e1 + 1):
                        genome_mapping[align.ref][pos] = 1
            for i in ns_by_chromosomes[align.ref]:
                genome_mapping[align.ref][i] = 0

    covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
    return covered_ref_bases, indels_info

This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.
This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).
From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast used_snps file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.


Software recommendations - DNA composition I'm looking to generate data based on the DNA composition of a region of my genomes (data is incomplete genomes from HiSeq runs in fasta format). I'm looking for software which will give me sliding windows for GC content, GC skew, codon bias etc. I would like to use something which I can use with or plug into R.

There are more R packages available that calculate GC content, for example Ape's GC.content() function.

For example:

> library(ape)
> data(woodmouse)
> GC.content(woodmouse)
[1] 0.3873347

With a sliding window is in Biostrings package.

> library(Biostrings)
> DNA <- DNAString("ACTGAAACCGTGGCAGTTTGAC")
> letterFrequencyInSlidingView(DNA, view.width=4,letters="CG")

For codon usage bias, I think cubfits or sscu packages might be useful.


how to export the monocle data into Seurat I have tried to export the monocle data into Seurat and got the error

"Error in Seurat::CreateSeuratObject(raw.data = data, normalization.method = "LogNormalize",  :
  unused arguments (raw.data = data, normalization.method = "LogNormalize", do.scale = TRUE, do.center = TRUE, is.expr = monocle_cds@lowerDetectionLimit)
Calls: exportCDS
In addition: Warning message:
In if (export_to == "Seurat") { :
  the condition has length > 1 and only the first element will be used"

Here is my code seurat_SAN_cds_subset <- exportCDS(SAN_cds_subset)

Is anyone familiar with this?

There's no default for the export_to option, so you just need to specify export_to="Seurat". Note that monocle2 doesn't seem to be compatible with the newest Seurat (3.1) releases.


How to show bootstrap values on a phylogenetic tree constructed with RAxML I apologize for this very basic question. With @Michael G.'s help, I have been learning how to use RAxML to build a phylogenetic tree (How to define an outgroup to build a robust amino acid tree) and after running my analysis I obtained several results files. I have both a RAxML_bestTree.result file and a RAxML_bootstrap.result file. How do I "combine" these files so that my tree shows the bootstrap values?

Any help will be greatly appreciated.

Depending on the search parameter settings RAxML produces several output files among which you normally find:

  RAxML_bipartitions.result: If you used the -f b option, this file
  will contain the input tree with confidence values from 0 to 100 drawn
  on its nodes! It is also printed when -f a x have been specified, at
  the end of the analysis the program will draw the BS support values on
  the best tree found during the ML search.

  RAxML_bipartitionsBranchLabels.result: Contains the same information
  as the file above, but support values are correctly displayed as
  Newick branch labels and not node labels! Support values always refer
  to branches/splits of trees and never to nodes of the tree. Note that,
  some tree viewers have problems displaying branch labels, they are
  however part of the standard Newick format.

from https://cme.h-its.org/exelixis/resource/download/NewManual.pdf

If you don't find the files in your output directory you did not specify the above parameters. If a rerun is not feasible you could you use the multinewick file RAxML_bootstrap.result with raxml or as input for the program splitstree to generate a consenus tree with different rules (majority rule, strict, etc.). However, these options don't guarantee to generate the topology of the maximum likelihood tree.


Phylogenetics analysis between diverse species Hope everyone is well.

So i was appointed with a task, and i need some guidance. Attention, the one that gave me the task may not have explained this very clearly. I will write what he exactly asked of me, and what he gave me as primary data.

So he gave me 4 bed files (human genome) with chr and coordinates(start/stop), and asked me to do: "sequencing/coordinates conservation between diverse species".

Is that doable with coordinates only or do i need sequences in fasta format (phast uses fasta format)? I just dont understand how i will find "conservation between diverse species" only with coordinates. It seems weird to me.

I have found these programs: TargetOrtho2, Phast and BUSCO(?). Are they any good? Should i continue with one of them that you may have experience with?

Thanks

I second commenter Devon Ryan's suggestion, I'm just going to give a little more detail about how to do that. You can plug the BED files into the UCSC browser following these instructions: https://genome.ucsc.edu/goldenPath/help/customTrack.html

Now you can look at the overlap between your bed tracks (black boxes I believe) and conservation (red box).

Assuming that is all you want to do, that should do it. If you want to do further data analysis or processing, you can download the conservation track as a .wib file from UCSC. You can convert this file to BED as described here, and use BEDTools to analyze it together with your BED of regions of interest.


ribo interval list for picard ensembl mm9 gtf Could some one please help me with understanding how to generate the ribosomal interval list that is required to when using picard metrics this step

java -jar picard.jar CollectRnaSeqMetrics \
      I=input.bam \
      O=output.RNA_Metrics \
      REF_FLAT=ref_flat.txt \
      STRAND=SECOND_READ_TRANSCRIPTION_STRAND \
      RIBOSOMAL_INTERVALS=ribosomal.interval_list

I have a couple of bam files that were aligned using mm9 ensembl gtf as annotation, i did a lot of search and found this link Ribosomal Intervals For Collectrnaseqmetrics but this link or any other source does not explain how to produce or find ribo interval list from ensembl annotation

At least with Ensembl annotations, you can grep rRNA foo.gtf > rRNA.gtf to get the annotated rRNAs. Note however, that the 45S cassette isn't in the mouse reference genome, so you'll never have alignments to 18S, 5.8S, or 28S rRNA. I always get the feeling that the Collect*Metrics programs are made with human sequences in mind, since the cassette is included in the human reference genome.

Anyway, while you can use the annotation and reference sequence to make picard happy, make sure to append the Rn45S sequence to your genome if you want to actually measure rRNA levels. In my experience this is only useful for rRNA-depleted samples, since the rates can vary wildly there.

Edit: To make a file compatible with picard's "interval list", one can do the following with the GTF file:

grep rRNA foo.gtf | awk 'BEGIN{OFS="\t"}{print $1, $4-1, $5}' > foo.bed

According to the picard documentation a BED file can be used as input and the above awk code extracts the appropriate columns and switches to 0-based coordinates used in the BED format.


Python equivalent to R GRanges and IRanges Problem:
I am trying to convert some codes written in R to Python and part of that conversion process is find classes equivalent to the GRanges and IRanges from the GenomicRanges R package in Python. https://bioconductor.org/packages/release/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesIntroduction.html#granges-genomic-ranges
I couldn't find any equivalent library written in Python that does the same operations as the aforementioned R classes. There is however pyranges library in Python but it is not as flexible as GRanges and doesn't have IRanges implemented in it. https://github.com/biocore-ntnu/pyranges
For example in R I could add the following in the GRanges class:
gr = GRanges(
seq = Rle(DF$chr),
ranges = IRanges(DF$Start,DF$End),
ref = DF$ref,
alt = DF$alt,
sampleID=DF$sampleID,
seqlengths = chromSize)

However, pyranges don't allow you to add arguments other than (chr, start, end, strand). However, in my case, I need to have the option to add sample_id, seqlengths and other arguments.
Thank you in advance.
Best

I am not sure if your question has anything to do with IRanges. If I get it correct, the limitation you have with pyranges now is adding extra columns to the meta data, like adding values to @values slot in granges.
You can do it using insert or setattr :
import pyranges
import pandas as pd
import numpy as np
gr = pyranges.random(10)

setattr(gr,'alt',np.random.choice(['A','T','G','C'],10))

gr.insert(pd.DataFrame({'ref':np.random.choice(['A','T','G','C'],10)}))

+--------------+-----------+-----------+--------------+------------+------------+
| Chromosome   | Start     | End       | Strand       | alt        | ref        |
| (category)   | (int32)   | (int32)   | (category)   | (object)   | (object)   |
|--------------+-----------+-----------+--------------+------------+------------|
| chr1         | 27431195  | 27431295  | -            | C          | G          |
| chr2         | 197045893 | 197045993 | -            | A          | G          |
| chr4         | 86316012  | 86316112  | -            | C          | A          |
| chr8         | 4560598   | 4560698   | -            | G          | G          |
| ...          | ...       | ...       | ...          | ...        | ...        |
| chr11        | 130395526 | 130395626 | -            | A          | A          |
| chr12        | 28557156  | 28557256  | +            | T          | G          |
| chr13        | 55519337  | 55519437  | -            | T          | C          |
| chr14        | 1807028   | 1807128   | +            | G          | A          |
+--------------+-----------+-----------+--------------+------------+------------+
Stranded PyRanges object has 10 rows and 6 columns from 10 chromosomes.
For printing, the PyRanges was sorted on Chromosome and Strand.

Or directly adding them like this (thanks to @ChrisRands for pointing it out):
gr.alt = np.random.choice(['A','T','G','C'],10)


TCGA data sets joining multiple files I have downloaded harmonised the TCGA-HNSC dataset and extracted it. It has 546 samples in total, so what I see in it are FPKM normalised values. The first column remains the same with 60000 rows and 2 columns that is first column ENSEMBL ID and second column FPKM value.

So now the issue is I want to join all the 546 txt file into one where the first column remains the same and rest 546 column would be joined. I'm using this command:

paste -d " " *txt > merged_file.txt

This is joining the all the data sets as ENSEMBL ID and FPKM value..

But I want 1st column as the ENSEMBLE ID and rest the sample FPKM any suggestion or help would be appreciated

In R:

flist = list.files("*.txt", path=".")
fl = lapply(flist, function(x) read.delim(x, header=F))
m = sapply(fl, function(x) x$V2)
row.names(m) = fl[[1]]$V1
colnames(m) = flist

That will give you an appropriately formatted matrix named m in R, which is probably where you want to continue processing things anyway.


Create GFF3 feature exporter. The information is below: Problem Statement
The GFF3 format is a commonly-used one in bioinformatics for representing sequence
annotation. You can find the specification here:
http://www.sequenceontology.org/gff3.shtml
The genome and annotation for Saccharomyces cerevisiae S288C is on the class server here:
/home/jorvis1/Saccharomyces_cerevisiae_S288C.annotation.gff

Note that this same file has both the annotation feature table and the FASTA sequence for the
molecules referenced. (See the '##FASTA' directive in the specification.)
Within the feature table another column of note is the 9 th , where we can store any key=value
pairs relevant to that row's feature such as ID, Ontology_term or Note.
Your task is to write a GFF3 feature exporter. A user should be able to run your script like this:
$ export_gff3_feature.py --source_gff=/path/to/some.gff3 --type=gene --attribute=ID --value=YAR003W

There are 4 arguments here that correspond to values in the GFF3 columns. In this case, your
script should read the path to a GFF3 file, find any gene (column 3) which has an ID=YAR003W
(column 9). When it finds this, it should use the coordinates for that feature (columns 4, 5 and 7)
and the FASTA sequence at the end of the document to return its FASTA sequence.
Your script should work regardless of the parameter values passed, warning the user if no
features were found that matched their query. (It should also check and warn if more than one
feature matches the query.)
The output should just be printed on STDOUT (no writing to a file is necessary.) It should have a
header which matches their query, like this:

gene:ID:YAR003W
.... sequence here ...

As an extra challenge, you can format the sequence portion of the FASTA output as 60-characters
per line, which follows the standard.
So far I have this below. I'm not even sure if this is correct to start with:
#!/usr/bin/env python3

for line in open("Saccharomyces_cerevisiae_S288C.annotation.gff"):
   line1 = line.rstrip()
   if line.startswith("#"):
       continue
   column = line1.split("\t")
   if len(column) != 9:
       continue
   id = column[8]
   type = column[2]

Can anyone help me with this? I'm totally lost. Thank you!

Your answer looks like a good start. You'll probably want to use sys.argv, argparse, or another argument parsing tool to read in options from the command line, then compare them to the fields (e.g. the "type" argument gets compared with the type column).
Given that attributes can be arranged in an arbitrary order, you'll probably need to read all of them from each line and store in a dict in order to determine whether the target is present. A regular expression search is another way to do that.
There are a few libraries available (e.g. from biopython) for reading and subsetting FASTA records.


Subset FASTA file by species name I have a problem: I've managed to download a massive fasta file of 1500 sequences, but now I want to split them into separate fasta files based on the genus.
EDIT
The fasta file looks like this:
terminase_large.fasta
>YP_009300697.1 terminase large subunit [Arthrobacter phage Mudcat]
MGLSNTATPLYYGQF...

>YP_009208724.1 hypothetical protein ADP65_00072 [Achromobacter phage phiAxp-3]
MSNVLLKQELDEWL...

>YP_009148449.1 large terminase subunit [Delftia phage RG-2014]
MSEPRKLVKKTLD...
...

I would like to end up with something like this:
mycobacterium_terminase.fasta
>Mycobacterium phage JAMaL
MVRKKPPPELE...

>Mycobacterium phage Bruin
MEVCGYTLDDI...

>Mycobacterium phage Zaka
MSLDNHLPELA...

salmonella_terminase.fasta
>Salmonella virus SETP7
MNVDITATEPQ...

>Salmonella virus SE2
MEGGSDSLIAM...
and so on...

So that the label on the alignment would have the name of the phage and not the protein's ID.
I've managed to extract the genus names of my organism with this:
outfile = open('species.txt', 'w')
with open('terminase_large.fasta') as fd:
    for line in fd:
        if line.startswith('>'):
            if '[' in line:
                name=line.split('[')[-1]
                name=name.split(' ', 1)[0]
                outfile.write(name[:] + "\n")
outfile.close()

And got to extract only the unique names with this:
lines_seen = set()
outfile = open('species2.txt', "w")
for line in open("species.txt", "r"):
    if line not in lines_seen:  # not a duplicate
        outfile.write(line)
        lines_seen.add(line)
outfile.close()

(Can I merge those two scripts together?)
Now, my genus names look like this:
Arthrobacter
Achromobacter
Delftia
....

I tried automating my script to get the Entrez data, but it gives me the 'Supplied id parameter is empty' message
My code looks like this:
from Bio import Entrez
Entrez.email = "example@example.org"

for line in open("species2.txt", "r"):
    searchterm = "(terminase large subunit AND viruses[Organism]) AND" +line+ "AND refseq[Filter]"
    searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=1000)
    searchResult = Entrez.read(searchResultHandle)
    ids = searchResult["IdList"]

    handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
    record = handle.read()

    out_handle = open('terminase_large_'+str(line[:-1])+'.fasta', 'w')
    out_handle.write(record.rstrip('\n'))

Can someone help me with it?

Splitting into multiple files and changing the IDs can be easily done:

perl -pe 'if(/>/){/\[(.*?)\]\s*$/; $_="> $1\n"}' file.fa |
    awk '(/^>/){name=$2} {print >> name".fa"}'

That assumes all your FASTA headers have [foo bar baz] as the last element of a line. It will create a file called foo.fa (the bacterium's name) with all sequences saved there.

Explanation

perl -pe : run the script given by -e on each line of the input file, and print the resulting line.
if(/>/) : if this line starts with a >. \
/\[(.*?)\]\s*$/ : match an opening bracket (\[), then capture (that's what the parentheses do, they capture a pattern so we can refer to it as $1) everything until the first ] (.*?\])
$_="> $1\n" : the $_ special variable in Perl is (in this case) the current line. So, $_=foo means "make the current line read foo. Since the -p prints each input line, changing the value of $_ means the changed value will be printed. So here, we are printing >, whatever was in the square brackets ($1) and a newline character.

The output of the perl command alone on your example input file is:

$ perl -pe 'if(/>/){/\[(.*)\]\s*$/; $_="> $1\n"}' file.fa
> Arthrobacter phage Mudcat
MGLSNTATPLYYGQF...

> Achromobacter phage phiAxp-3
MSNVLLKQELDEWL...

> Delftia phage RG-2014
MSEPRKLVKKTLD...

So, we now pass it through awk which does:

(/^>/){a=$2} : if this line starts with an >, save the second field (the bacterial species) as the variable name.
{print >> name".fa"} : print each line into a file whose name is the current value of the variable name with a .fa. extension.

If you prefer python scripts to the one-liner approach, you can do the same thing with:

#!/usr/bin/env python
import re
outFile = None
with open("file.fa", "r") as inFile, open(species + ".fa", 'a') as outFile:
    for line in inFile:
        line = line.rstrip()
        if line.startswith('>'):
            regex = re.compile('.*\[((.+?)\s+.*?)\].*')
            matches = regex.search(line)
            species = matches[2]
            outFile.write('>%s\n' % matches[1])
        else:
            outFile.write("%s\n" % line)

As for your script, you've got the right idea, but have a small bug. You forgot to remove the \n from your input file, so it looks for Arthrobacter\n instead of Arthrobacter. The golden rule of debugging is "print all the things". If you add print("Searchterm: ",searchterm)  you will see:

Searchterm:  (terminase large subunit AND viruses[Organism]) ANDArthrobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDAchromobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDDelftia
AND refseq[Filter]

So, you need to remove the newline characters and add a space like so (I also made it a bit more "pythonic" and conforming to the Python syntax guidelines):

#!/usr/bin/env python
from Bio import Entrez
Entrez.email = "example@example.org"

with open("species2.txt", "r") as in_handle:
    for line in in_handle:
        line = line.rstrip()
        searchterm = ("(terminase large subunit AND viruses[Organism]) " +
                      "AND %s AND refseq[Filter]" % line)
        print("Searchterm: ", searchterm)
        searchResultHandle = Entrez.esearch(db="protein",
                                            term=searchterm, retmax=1000)
        searchResult = Entrez.read(searchResultHandle)
        ids = searchResult["IdList"]

        handle = Entrez.efetch(db="protein", id=ids,
                               rettype="fasta", retmode="text")
        record = handle.read()

        with open('terminase_large_' + str(line[:-1]) + '.fasta', 'w') as out_handle:
            out_handle.write(record.rstrip('\n'))


Retrieving gene matrix from internet in GSEA desktop This error appeared to me, and I have the program connected to the internet. I already installed java 8, but it still doesn't get me the gene set databases.

Have you checked firewall settings on your laptop?

See Firewall / FTP connection issues at:

GSEA wiki - known issues


WGCNA module preservation analysis I have case control studies,4 different condition where I have 1 control and against which I'm have 3 comparisons to make.

I did WGCNA for individuals samples .Here I have the kept the number of genes equal in all the cases .

Now my question is whether the number of modules which are obtained from WGCNA would be similar in number if I run module preservation between two condition .

Lets say i get 15 modules in Control condition and 16 treated condition if I run module preservation the number of modules which would come would be 15+16= 31 or its independent .

The power I have kept same for both WGCNA and module preservation analysis. But Still I get different module numbers.

Any suggestion would be highly appreciated

Your question is not clear, you seem to mix the relevant terminology...

First, you cannot run WGCNA for individual samples, only for individual conditions (with or without WT, depending on what question you're trying to answer and assuming you have enough samples in each condition). If that's what you have done, fine. The number of modules in each independent analysis will generally be different, that's because each condition is different. The module identification built into WGCNA as the default approach does not allow the user to specify the number of modules beforehand, but if you insist on getting a pre-specified number of modules, you could run WGCNA using the "manual" way (for example, Tutorial I, section 2b) and substitute another clustering method (e.g., PAM) that allows you to specify the number of clusters. You may find the results inferior though.

Second, module preservation does not result in new modules. Module preservation evaluates whether modules found in data set A are preserved in a different data set B. To get new, common modules, you could run a consensus network analysis on pairs of conditions or all 3 conditions together. The number of modules you would find in a consensus analysis cannot be predicted from the number of modules you find in individual data sets.


Import gene list to Seurat to define cell types Is there a way to import gene list into Seurat to define cell type? The default cell types in Seurat is not enough for our research. For example, we want to mark a subtype of B cells in Seurat, but seems like Seurat only have general B cell, I'm planning to download gene list from Genomic Cytometry, and then import the gene list into Seurat and then Seurat define cell types in the cluster based on the gene list I import. Does Seurat support this function?

Seurat does not define cell types by name. It clusters and assigns each cell to a cluster, from 0 to X.
If your data has the cell type (e.g. B,T, Mast cells) it means that someone annotate the clusters so that they have a biological meaning.

You can assign different names to the clusters by using the AddMetaData function.

sampleID  cellType
cell1     B
cell2     B
cell3     T
cell4     T

newMeta <- oldMEta
newMeta$origType <- newMeta$cellType
newMeta$cellType <- c("B1","B","Treg","Tkill")

object <- AddMetaData(object, newMeta)

The new metadata should look like

sampleID  origType  cellType
cell1     B         B1
cell2     B         B
cell3     T         Treg
cell4     T         Tkill

You can then swap the labels using SetAllIdent(object = test, id = "cellType"). Your tSNE or uMAp, will now show the new cell types.

Have a look at Manually define clusters in Seurat and determine marker genes


Interpreting bayesian phylogenetic model using Tracer? I am doing phylodynamic analysis of virus sequences using bayesian method using BEAUti and BEAST. As a part of my analysis I am using Tracer to check how my model is working.

I can interpret some of the results like from the Tracer like :

clock.rate - rate of nucleotide substitution

treemodel.rootheight- time when divergence occured

ESS- measures sampling quality (ideally has to be greater that 200)

But there are other many terms on this analysis which I have no idea how to explain.

Please can anybody help me to find out what other aspects of the evolution can be explained using this table.

posterior   -33531.076
prior   -481.982
likelihood  -33049.094
treeModel.rootHeight    32.271
tmrca(Brazil)   3.853
tmrca(Canada)   16.644
tmrca(China )   11.752
tmrca(Columbia) 1
tmrca(Thailand) 1.093
tmrca(USA)  32.271
tmrca(allseq)   32.271
constant.popSize    43.595
CP1.kappa   18.716
CP2.kappa   17.16
CP3.kappa   22.796
frequencies1    0.224
frequencies3    0.221
frequencies4    0.26    5001
CP1.mu  0.76
CP2.mu  1.221
CP3.mu  1.019
clock.rate  1.354E-3
meanRate    1.354E-3
CP1.treeLikelihood  -9462.258
CP2.treeLikelihood  -12520.305
branchRates 0E0
coalescent  -467.376

When I have a moment I will fill in the blanks

You must have an RNA virus (nothing else will mutate that fast). The basal group and root is the USA. It estimates the time to that basal event is 32.3 (likely years).

"frequencies" frequency of ACGT where the missing one will be 0.25.

Its a simple model of constant population size.

ESS is the most important parameter to ensure convergence of the Bayesian simulation .

Likelihood is the second most important parameter. I'll edit again shortly.

I'll fill in one more blank and summarise ...

Kappa ...
This will be the transition transversion ratio in HKY. I suspect it also will be transition transversion ratio in Kimuras 3-parameter model, this is because ... the transversions are further split into two mutation groups, ie. purine to pyrimidine has 2 further combinations and you end up with 2 ratios.  that makes sense, transition/transversion (kappa 1), transition/transversion 1 (kappa 2), transition/transversion2 (kappa 3) ....

CP1.kappa 18.716
CP2.kappa 17.16
CP3.kappa 22.796

CP2 and CP3 are either side of CP1 .. which is what you would expect on my interpretation.

These are very high values by the way, essentially there are very few transversions in this data, so ... I dunno. You need to check the data, if you see lots (and lots) of A-G and T-C mutations, then thats what I believe it is saying.

Briefly posterior and prior are the fundamental method of a Bayesian calculations, the output is probability expressed as a log, so the more negative the higher the probability (the prior looks like a random tree, or a reasonable tree with really rubbish parameters, likely the former).

Summary
Essentially I am going to call it a day here. The real issue is that the output summarise the entire modern theory in phylogenetics, so its a very long response required to address all of the parameters and I honestly don't think it helps your analysis. What I am clear about is:

* You've omitted a gamma distribution (I think Tracer would summarise this) *

So my advice is:

ESS MUST be homogeneous between runs - big trouble if not
gamma distribution MUST be run

Beast can be tricky getting gamma to converge via ESS, just a word of warning.

I can't see the full robustness, because its absent... However, from the output, I think its pretty clear, that you (well me) can more or less see the model you ran.


Why is FPKM still used for gene expression studies? Perhaps my understanding is wildly misguided, but I'm seriously confused about why people are still interested using FPKM values for cross sample gene expression analysis. My understanding is, that for some gene ${ i }$, its $FPKM$ is calculated as follows:

$$ FPKM_i \propto \frac{ \mathrm{fragments\ mapping \ to\ gene \ }i } { \sum_j \mathrm{ fragments\ mapping\ to \ gene\  } j } \bigg/
\mathrm{effective\_length}(i) $$

For comparing gene expression within a single sample, this doesn't seem too harmful to me. In fact, this seems like a good idea.

However, the use of $FPKM$ or even the famous $log_2(FPKM+1)$ in a similarity analysis seems to me plain wrong. True differences in RNA abundances of a couple of genes would significantly impact the FPKM values of the other genes leaving  perhaps confused about the true differences between the samples.

Am I missing something here? Is it ever meaningful to use FPKM values for analyzing across samples?

Is it ever meaningful to use FPKM values for analyzing across samples?

One should never use FPKMs for anything important. They can occasionally be useful for plotting, but even in that case one needs to construct the FPKMs from properly normalized data and not use the original definition.

It's becoming ever rarer to see people using FPKMs these days, thankfully.