text:
Aligned base strand in pysam `pileupcolumn`
Is it possible that on the same pileupcolumn object (i.e a specific aligned base) I find the same base but aligned to different strands? For example, can I find a coverage of 50 'A's that align to strand '+', and on the same time 50 'A's that align to strand '-' within the same index? I expect that if 'A' is aligned to strand '+', I will find 'T' aligned to strand '-'.

Please [edit] your question and clarify what you are doing. What code are you running? On what input? Anything is possible, including mismatched bases either in the real data or because of a sequencing error, but we can't know if you don't give us more details.

As @terdon mentions, it's certainly possible to have an A aligned on both strands. If a read has a mismatch at the location of interest and that mismatch is the complement of the reference base, then you are likely to see A in both strands. However, it is very unlikely you will see the same base at high coverage in both strands. For example, if your total coverage is 55, and you see 50 As on the + strand and 50 As on the - strand, it's likely that the pileup object is reporting the reverse complement sequence of the reads that align to the - strand for your convenience.


text:
Where to assign a read pair when a mate spans a HiC bin boundary?
Chromosome capture experiments, like Hi-C, typically bin the genome into distinct bins. A contact matrix, $M$, is given by $M_{i,j} = |\{pairs | mate_1 \in bin_i, mate_2 \in bin_j\}|$.
If one of the read mates spans a bin boundary, where should that pair be assigned?
For example, if $mate_1$ overlaps $bin_i$ and $mate_2$ overlaps both $bin_j$ and $bin_{j+1}$, should this read pair be assigned to $M_{i,j}$ or $M_{i,j+1}$? Should it be randomly assigned to either bin with equal probability? Should it be assigned to both?
I know this may be a rare event, but if Hi-C datasets are generated using billions of reads (such as Rao et al., Cell, 2014), this is bound to happen.

[Old question, but possibly still relevant to answer.]
This will not necessarily be rare. One or the other of the mates could span the ligation junction itself, in which case your read start and end could be entirely different chromosomes. In fact, some tools specifically measure the quality of Hi-C libraries by measuring the number of such ligations directly observable within reads. This will be a much larger number of reads than bin-crossing reads.
The safest thing to do in all such cases (IMO) is to take the read start position as the location of interest for each mate. That way, if a junction does occur, you still have an unambiguous position which is your best guess of where each end of the ligation product actually originates. This gets around both the minor issue of the bin-crossing problem and the major issue of the ligation junction problem. There is probably some imprecision here and you throw out some information, but sorting out the complex possible situations otherwise seems like a lot of hassle for minimal gain.

I think what I was referring to here were undigested contiguous reads, like Figure 1 of the HiCUP method paper describes. In this case, HiCUP marks this as an invalid read pair and filters them out. From my experience, there weren't enough of these occurrences (thankfully) that including or excluding them would have affected the contact matrix that much. At the time I asked this question, the qc3C method you referenced wasn't published yet, but I think that's a good reference to link here, now that it is available.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
how do I convert a vcf file for a GWAS study (using R package vcfR)
I have a vcf file with individual level genetic data, that I read using the R package vcfR
gen <- read.vcfR('test.vcf')

Scanning file to determine attributes.
File attributes:
  meta lines: 9
  header_line: 10
  variant count: 152
  column count: 3798
Meta line 9 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 152
  Character matrix gt cols: 3798
  skip: 0
  nrows: 153
  row_num: 0
Processed variant: 153
All variants processed

How do I convert this file to a data frame that I can use for a GWAS or Mendelian Randomization study?
I was told that this should be a data frame where each column is a genetic variant and it take values 0, 1, or 2, depending on the number of "mutations".
I never worked with genetic data before and do not know much about genetics.

Welcome to Bioinformatics.StackExchange! GWAS is a broad area, and there are numerous tools you could use. The format of your data will depend on which tool you are using, and hopefully this will be clearly documented in the tool's manual. If on the other hand you're trying to implement your own tool, you'll have to be much more specific about what you want (or are required) to accomplish.

Welcome to the site! Could you please explain what have you tried? For instance, what happens if you use as(gen, "data.frame")?

I found that I can use extract.gt, then transform the matrix. Entries in this matrix are "0/0", "0/1", "1/0", and "1/1". I now have to recode these to 0, 1, 1, 2.

If you want to do GWAS or mendelian randomisation, you can do it with Plink (v2) which should be faster than R.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How do I get the adjacency matrix from gene similarity network?
I am exploring application of Graph convolution network on gene similarity network. I am trying to explore
genemania.org and https://string-db.org/
Please let me know where can I find the data for my research in gene interaction network, which allows me to download the data in the form of an adjacency matrix.

Welcome to Bioinformatics.StackExchange! To help us help you, it's always a good idea to let us know where you've looked and/or what you've tried so far.

Also, please [edit] your question and show us an example of your desired format. Do you only want the names of interacting proteins? Do you want some sort of value for the edges? Do you want all interactions (STRING has a lot of noise) or only high quality one? Maybe you only want binary, direct protein-protein interactions (so no immunoprecipitation results, or no complexes for example)? The more details you give us, the better your chances of getting a useful answer.

The STRING download page is here. I am not sure that they provide any data in the form of an adjacency matrix, as the number of nodes (rows/columns) would be quite large, and for this reason they are likely to represent the data as lists of links, for compactness.
Their data will likely have to be cleaned for any purpose you are interested in. You will have to make some decisions on how to filter the data. For example, as suggested in the comments, you may want to filter on the quality of the evidence supporting interactions, or based on species of the protein, or other things.
It's possible that you may have better luck going with another database that has done some of these steps already. For example, there are relatively compact human interactome datasets here, or yeast interactome datasets here. In this case you will still need to convert from a list of links to an adjacency matrix, however this is not very difficult if you can use a scripting language.
Notably, you may have to figure out how to distinguish missing data from tested non-interactions in these cases (or decide that missing data is ok). In most cases people will not have exhaustively tested every cell of the all-by-all interactome adjacency matrix.
Hope that helps.


text:
Why is it necessary to add hydrogen and delete water before protein-ligand docking?
What is the reason for adding hydrogen and removing unnecessary water molecules from the protein structure before protein-ligand docking? FYI, the tools I used for docking is GOLD.

I think if you could add some more context it would be helpful to answer. Who said it is needed? Under which conditions are you looking for a docking? Sometimes  water is need in order to enable interaction or form stable bounds (AFAIK)

No one say that but from the software, it say "delete water IF REQUIRED". I just don't understand what it meant by "IF REQUIRED" and when it is necessary to remove water or add hydrogen. Btw, from your reply, does it meant that sometime, water is required for an enzyme to function?

Yes, sometimes water is needed for an enzyme to function. When it is necessary to remove water is a hard question that depends on what environment are you simulating and the protein you are modelling

My answer is mainly based on molecular dynamics,but in this context it shouldn't make much of a difference.
The main reason has to do with the force-field-water model used during the simulation. Examples of water models are SPC, TIP3P or TIP4P; same molecule (water) but different topology.
A pdb usually comes from different sources and it can contain molecules that are not standard and therefore not recognised by the model of your simulation. Usually, once the structure is loaded, your system is re-built so that it conforms to the model topology. There are exception though, for example if your water is part of a functional active site. In that case you can keep it, but its structures (topology of the water molecule) has to be recognised by the model.
Hydrogens are usually added to "saturate" and stabilise the bonds.

Totally correct for MD. However, a wee note: most ligand docking program use implicit water (except special case such as WaterDock), which means that TIP3P waters or similar are not added, but a fudge factor is used to account for the dynamic but stabilising polar properties of water.

Hydrogens
Most crystallographic structures are X-ray and the hydrogens are absent. In the case of neutron diffraction models, there may be many missing high b-factor ones (example in PDB:5MOP). As a consequence models are reprotonated. Also missing carbons and missing loop residues are added for some in silico applications too (not docking).
However, it is always good to manually tweak the protonation. In the case of enzymes with protanatable side chains in the active site, some reading is required to get the charges right: by default aspartates will be unprotonated, but they may be protonated (glutamic acid) in the active site and it is important to get it right (cf. this post about it). Some times its the ligand parameterisation may give the physiologically incorrect protonation (cf. this post).
Waters
Crystallographic waters can be important catalytically and structurally, but can be pushed out and actually it is a good thing: you have a hit that take its place. Polar bonds are stronger than hydrophobic ones and most types of π bonding.
Waters in most docking software is implicit and is modelled as a parameter in the force-fields —think of it as a ghostly watery haze. So a left behind water, not only will get in the way of the docked small molecule, it may for some programs cause issues as they will try to dock it. A notable exception is WaterDock mod for AutoDock, but the jury is out on whether it is better. Water is polarisable in reality in a very complex way (protonated (OH3+) or deprotonated (OH-) with very different charges) and most water models do not cut it and instead do not account for this.
The way implicit waters are computed differs in some applications, which is also something to look into.

The molecular docking programs usually tread the target structure as rigid. It does not know which atom belongs to the receptor and which is bulk water. So, if you keep the water in the structure the program cannot properly place the ligand there.


text:
scRNA-seq multi-dataset integration for small datasets
There have been a few methods proposed for integration (or batch correction) of scRNA-seq datasets, such as Seurat CCA, MNN Correct, Scanorama, and Harmony. The concern is generally about the maximum number of cells that they handle, but I haven't seen any discussion about the minimum number of cells. I am confident they can all handle 10k cells reasonably well and will fail with 10 cells, but where do you draw the line? Is there a method that works best for small datasets?
For example, with plate-based platforms like Fluidigm, many experiments only have 96 cells and potentially much less after quality filtering. How can those be used?

I hope my answer is sufficient speculation as not many papers have been written specifically on this topic.

TLDR: Harmony can work on 106~ cell samples but has less versatility then methods like BATMAN. BATMAN is useful if you need your data for differential gene expression and single-cell eQTL and can work on ~200 cell data.
Well their certainly is a minimum which comes down to several factors. As many of these rely on algorithms like KNN and other clustering algorithms for integration the threshold is dependent on the quality of the data and the  algorithms abilities to predict reasonable clusters from that data. As well as the downstream uses of the data...
BATMAN highest performance on simulated set 200 cells and 1000 cells
With 200 cells from one batch and and 1000 in the other only BATMAN succeeded in clustering the high dimensionality simulated data well. If you need your data for differential gene expression and single-cell eQTL, BATMAN seems to be the correct choice, for small sets. "The downside of these methods is that they operate in latent space, which limits their interpretability and use in downstream analyses such as differential gene expression and single-cell eQTL analyses."
"A shows that the two original datasets become better mixed with each other. BATMAN has the highest performance not only when considering the top two principal components but also when considering more top principal components. It is the only method that manages to efficiently maintain a high iLISI score in a larger number of dimensions."
Harmony is best if... differential gene expression and single-cell eQTL is not necessary
"Moreover, we show that Harmony requires dramatically fewer computational resources. It is the only available algorithm that makes the integration of ~106 cells feasible on a personal computer."
Computing time benchmarks-section
It seems as the runtime of these algorithms is less than 2 minutes, using a set with less than ~2000 cells will not benefit your computational time much.
"When processing a small data set of ~ 2000 cells, all four methods took less than 2 min."
"To obtain such datasets, we downsampled the MCA and TM datasets to obtain a total of 9 sets of data containing between ~ 2000 and ~ 140,000 cells, while the number of highly variable genes (HVGs) was controlled in a range from ~ 2000 to ~ 3000 (Table S1)."


text:
Help with 1D^2 library shearing
I've done my 1st trial of 1D^2 whole genome sequencing using LSK-SQK308 kit with R9.5 flowcell. Without doing library fragmentation, I encountered the unexpected library shearing during sequencing.

But I didn't see significant portion of low-molecular-weight fragments in the QC gel photo before and after library preparation.

Anyone using nanopore has encountered this before? Any explanation? And how to avoid it? I could not find evidence of the problem in my QC steps...
Thanks in advance!

Short template DNA is a big problem with nanopore sequencing; it doesn't easily show up on a gel if it's well distributed across different lengths (as appears to be the case here from the read length plot), and sucks up both adapter and sequencing capacity.
The usual way to deal with it (for reads in the size ranges seen here) is to reduce the bead concentration when doing magnetic bead cleanup. Try something like 0.6X beads (e.g. 30μl of Ampure XP beads added to 50μl template DNA solution), which should still allow most of the long template sequences to bind, but leave shorter reads in solution.


text:
How to predict stop codons in Illumina reads?
I have Illumina MiSeq paired-end reads from 150bp amplicons mapped to my reference genome (> 1000X coverage). These reads have indels that may or may not induce frameshifts. If the indel induces a frameshift (i.e. indel is not a multiple of 3bp), I want to know if that leads to a stop codon later in the read. Is there any package/method that would allow me to do this?
Few extra details

I might be wrong but I think this is different than gene prediction as I would want the tool to scan each read in the correct frame; not try different frames to detect ORFs etc.
I expect mosaicism, so I do not want to compute a single assembly/consensus and predict stop codons in it.

Here is how I imagine such a script/tool would work:

Take the position of the first translated codon of the exon I am interested in (from Ensembl for instance?). Eg. gene abc/exon1 starts at chr1:1000.

Then for each read:

Look at which position the read is mapped to my reference. Eg. read 1 is mapped to chr1:1050.
Compute the position of the first full codon in the correct frame in my read. Eg. start codon at chr1:1000 but my read starts at chr1:1050, the first full translated codon in my read should start at chr1:1052 (I think?).
Split all my read into 3-mers starting at that position and build something like a 3-mer/position dictionnary.
Look for stop codons and report the positions.

Does it exist?

Could you give us a bit more context? Is this targeted sequencing of a eukaryote? Or whole genome bacterial sequencing? So as input you have a reference genome (FASTA), annotated/predicted genes (GTF/GFF3), and your NGS reads (FASTQ/BAM)? And you'd like to predict stop codons based on the read mappings? Do you have reason to believe the annotated stop codons are incorrect? Are you looking for nonsense mutations?

@DanielStandage It is targeted sequencing of small (150bp) amplicons from zebrafish genomic DNA. I am sequencing after inducing mutations, i.e. my reads contain various indels. This is why I need to do this for each read in my alignment. I am looking for nonsense mutations, but not only. Eg. there is a 2bp deletion early in the read, I want to know if there will be a stop codon downstream as a result when the RNA gets translated.

Let us continue this discussion in chat.

I just recently discovered the ampliCan R package from the Valen lab after asking this question. I highly recommend it. It does all of what I was expecting & more.
To install: https://bioconductor.org/packages/release/bioc/html/amplican.html
To run (& other explanations): https://rdrr.io/bioc/amplican/f/vignettes/amplicanOverview.Rmd
This is the original paper & Do not forget to credit the people!
https://genome.cshlp.org/content/early/2019/03/08/gr.244293.118.abstract

With 150bp reads, you might be able to blastx all the unique sequences you've got against the target proteins.


text:
load the phenotype data for ballgown
I am trying to reproduce the work of this paper
[1], and I have run StringTie successfully, but after that I have to run Ballgown but could not understand this command:
>pheno_data = read.csv("geuvadis_phenodata.csv")

Actually, I dont know how to make the .csv file of my samples. What phenotype information should this csv contain?
I have checked the supplementary data of this research article but didn't find any .csv file.
How I can make the .csv file for my samples to load the phenotype data for ballgown?

Load the phenotype data for the samples. An example file called
  geuvadis_phenodata.csv is included with the data files for this protocol
  (ChrX_data). In general, you will have to create this file yourself. It contains
  information about your RNA sequencing samples, formatted as illustrated in
  this csv (comma-separated values) file. Each sample should be described on
  one row of the file and each column should contain one variable. To read this
  file into R we use the command read.csv. In this file the values are
  separated by commas, but if the file were tabdelimited you could use the
  function read.table.

pheno_data = read.csv("geuvadis_phenodata.csv")

Pertea et al, Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, Nature Protocols volume 11, pages 1650–1667 (2016)

You can easily do it via Notepad++ or even Excel. Just don't forget to save the output file as a .csv file (comma-separated value).
The info inside should look like (my working example):
id,state
lar_control141,control
lar_control142,control
lar_control143,control
lar_probiotic142,probiotic
lar_probiotic143,probiotic
lar_probiotic144,probiotic
lar_probiotic145,probiotic
lar_probiotic146,probiotic

That should work


text:
Extracting expression data in r from GEO and got a S4 object
I want to open the expression data GSE9838 dataset from GEO using GEOquery R package.
gset <- getGEO(filename="GSE9838_family.soft.gz")

The result of this code said it's S4 object of class GSE.
How can I get the expression data of GSE9838?

Welcome to the site. Have you looked for the tutorial of GEOquery? You can see them with vignette(package="GEOquery") . Also have you looked around on this site and in support.bioconductor.org for how have other people done this?

Possible duplicate of Extracting expression data from GSE dataset downloaded from GEO

Not sure if that method works in this case

@Oka Why shouldn't it work? Have Qiran tested it? At least mentioning that it didn't work it would help the people trying to answer the question (and future readers)

@llrs I´ve checked if it would, and it didn´t (at least on the first trial). The method worked fine with GSE used in the example, but didn´t seem to work in this case for some reason..

Ok, how it didn't work, could you post the code so I could see what is the error message and the options you used?

datExpr = exprs (data[[1]]) will just give an object with one header line (      GSM248431 GSM248432 GSM248433 GSM248434 GSM248435 GSM248436 GSM248437 GSM248438 GSM248439 GSM248440 GSM248441 GSM248442)

Well, if you use a SOFT formatted files:
#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that "Raw data provided as supplementary file. Processed data is available on Series record". So one may want to check series record instead:
# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz)

# and check what´s inside:
exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?
If you want to download them, you can:
filePaths = getGEOSuppFiles("GSE9838")
filePaths
untar("GSE9838/GSE9838_RAW.tar", exdir = "GSE9838/CEL")
list.files("GSE9838/CEL")

...what will allow you to see all the CEL files with the data.

Thank you so much! I will try to open the raw data.


text:
What is the difference between a Bioinformatics pipeline and workflow?
I want to understand the difference between pipeline systems and workflow engines.
After reading A Review of Scalable Bioinformatics Pipelines I had a good overview of current bioinformatics pipelines. After some further research I found that there is collection of highly capable workflow engines. My question is then based on what I saw for argo. I would say I can be used as a bioinformatics pipeline as well.
So how do bioinformatics pipelines differ from workflow engines?

Great question! Note that from a prescriptive standpoint, the terms pipeline and workflow don't have any strict or precise definitions. But it's still useful to take a descriptive standpoint and discuss how the terms are commonly used in the bioinformatics community.
But before talking about pipelines and workflows, it's helpful to talk about programs and scripts. A program or script typically implements a single data analysis task (or set of related tasks). Some examples include the following.

FastQC, a program that checks NGS reads for common quality issues
Trimmomatic, a program for cleaning NGS reads
salmon, a program for estimating transcript abundance from NGS reads
a custom R script that uses DESeq2 to perform differential expression analysis

A pipeline or a workflow refers to a particular kind of program or script that is intended primarily to combine other independent programs or scripts. For example, I might want to write an RNA-seq workflow that executes Trimmomatic, FastQC, salmon, and the R script using a single command. This is particularly useful if I have to run the same command many times, or if the commands take a long time to run. It's very inconvenient when you have to babysit your computer and wait for step 3 to finish so that you can launch step 4!
So when does a program become a pipeline? Honestly, there are no strict rules. In some cases it's clear: the 10-line Python script I wrote to split Fasta files is definitely NOT a pipeline, but the 200-line Python script I wrote that does nothing but invoke 6 other bioinformatics programs definitely IS a pipeline. There are a lot of tools that fall in the middle: they may require running multiple steps in a certain order, or implement their own processing but also delegate processing to other tools. Usually nobody worries too much about whether it's "correct" to call a particular tool a pipeline.
Finally, a workflow engine is the software used to actually execute your pipeline/workflow. As mentioned above, general-purpose scripting languages like Bash, Python, or Perl can be used to implement workflows. But there are other languages that are designed specifically for managing workflows. Perhaps the earliest and most popular of these is GNU Make, which was originally intended to help engineers coordinate software compilation but can be used for just about any workflow. More recently there has been a proliferation of tools intended to replace GNU Make for numerous languages in a variety of contexts. The most popular in bioinformatics seems to be Snakemake, which provides a nice balance of simplicity (through shell commands), flexibility (through configuration), and power-user support (through Python scripting). Build scripts written for these tools (i.e., a Makefile or Snakefile) are often called pipelines or workflows, and the workflow engine is the software that executes the workflow.
The workflow engines you listed above (such as argo) can certainly be used to coordinate bioinformatics workflows. Honestly though, these are aimed more at the broader tech industry: they involve not just workflow execution but also hardware and infrastructure coordination, and would require a level of engineering expertise/support not commonly available in a bioinformatics setting. This could change, however, as bioinformatics becomes more of a "big data" endeavor.
As a final note, I'll mention a few more relevant technologies that I wasn't able to fit above.

Docker: managing a consistent software environment across multiple (potentially dozens or hundreds) of computers; Singularity is Docker's less popular step-sister
Common Workflow Language (CWL): a generic language for declaring how each step of a workflow is executed, what inputs it needs, what outputs it creates, and approximately what resources (RAM, storage, CPU threads, etc.) are required to run it; designed to write workflows that can be run on a variety of workflow engines
Dockstore: a registry of bioinformatics workflows (heavy emphasis on genomics) that includes a Docker container and a CWL specification for each workflow
toil: a production-grade workflow engine used primarily for bioinformatics workflows

Thanks, Daniel. Your answer outlines very well the details about the bioinformatics workflows. Your argument, that a very high level of software engineering is needed to use my listed workflow engines, is on point . If I understand correctly generally speaking bioinformatics pipelines can be pressed into workflow engines. So thanks again. That helped a lot.

Nice answer! I'm biased, but I feel like nextflow is also worth a mention :) And also WDL to complete the picture of common workflow tools. There's a nice "awesome-pipeline" repo list of workflow tools here too: https://github.com/pditommaso/awesome-pipeline - a little more comprehensive than the one in the original post, perhaps with more of a bioinformatics slant.

Thanks for the comment @tallphil. I too am biased by the tools with which I have experience. :-) Feel free to edit my answer to add a bit about nextflow and WDL.


text:
Significant gene set testing - limma
While I have used my title, it is more of a general statistics question regarding how these significance tests are carried out, these include decideTests and mroast / roast among others.
To give a brief overview of my experimental design, I am analysing bulk RNASeq datasets and there are two main factors in the design matrix: patientID and casetimeViremia; where the latter identifies whether a patient is viremia positive or negative at a particular time point (1-4). This leaves us with 5 different possibilities for casetimeViremia: 1-min, 2-plus, 3-min, 3-plus, 4-min.
The design matrix was created in R by running:
model.matrix(~patientID + casetimeViremia, data = dgeObject$samples)
Where dgeObject is an object of the function DGEList.
For the purpose of answering this biological question, we wish to conduct the differential expression testing by seeing whether all the previously mentioned casetimeViremia timepoints are separately and significantly expressed relative to each person's first timepoint, which is labelled as: 1-min.
Now, I understood that should I wish to conduct this test, I do not need to specify a contrast matrix in the mroast function and can just do this by calling the column number of interest.. But in that case, how does mroast know that I wish to do the contrast relative to the 1-min timepoint? And, just wondering, should I wish to define a contrast matrix, how would I do this given that there are two colnames missing from the design, in this case one patientID and the 1-min timepoint?
Thank you all for your time.

The mroast function has an argument to specify which contrast do you want to test, quoting from the help page:

contrast
  contrast for which the test is required. Can be an integer specifying a column of design, or the name of a column of design, or a numeric contrast vector of length equal to the number of columns of design.

So if you have the contrast of 1-min timepoint positive vs 1-min timepoint negative you can use this.
Read the help pages in detail, specially the help pages in limma are quite extensive and detailed to answer most questions if not all about how to use the software.

thanks for your response! I have read the help pages and they were of no help in understanding this process. I even tried looking at the source code for the function but it was not available.

what exactly do you mean by "if you have the contrast of 1-min timepoint positive vs 1-min timepoint negative you can use this."? I understand that I can call the column 2-plus however, how does just calling on this column within the function specify that the contrast must be done relative to the 1-min timepoint?

A contrast is a comparison if you use a colum of the design you can compare 1-min timepoint vs the others but I assume that this is not what you want. It would be helpful if you posted your design matrix and your question (and your differential analysis )

So are you saying that should I call on any of the other columns in the design, it would default the contrast to {thatTimepoint - 1-min} ? This is in fact what I want, I just want to understand the mechanism behind it. And, hypothetically speaking, how I could run other contrasts should I wish to in the future. Also, good point, I will upload part of the design matrix in a few hours to give a better idea. Thanks again.

If you want to understand the mechanism look the code. I didn't mention the default of the contrast, and in fact it is the last column of the design matrix

Could you upload a link the source code if you know where it is? I have looked all over and could not find it.

It is in your computer but you can also access it from the tar.gz file from bioconductor page.

Yes I am familiar with these files.. Unfortunately though it does not contain the mroast() function in them.. Thanks anyway though for your help!

The mroast function is in the R/geneset-roast.R file (I found via grep -n "mroast" *.R. If this answer has helped you, please consider accepting this answer. If not please improve the question as I explained some comments above.


text:
Displaying soft-clipped nucleotides in samtools tview
There are nicer genomics visualization tools available, but the samtools tview command is almost always my go-to for a quick first look at read alignments. I just brought up the following locus in tview.
81055951  81055961  81055971  81055981                                                81055991  81056001  81056011  81056021
GTGGAGGTCGGCAGCAGAGCCTGTGGCAGCT**********************************************CTGAGGGTCGGCGTGGCCTCCTGGTGGAAGTTGCACCTGGTTTGCTT
...............................                                              ...............................................
...............................**********************************************..       ......................................
   ,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,       .....T................................
.....                  ,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
.......                      ,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,                                                                    ...............................................
.............                                                                .........................................A.....
........................                                                     ...............................................
...............................**********************************************..         ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,           ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************.A......                 ,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************..                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,    ,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,        ,,,,,
...............................**********************************************..                               ,,,,,,,,,,,,,,
...............................**********************************************..
...............................**********************************************..
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..................T............**********************************************...............................................
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..............G................**********************************************...............................................
,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
...............................CTAAGCAAGTAACTGTGGTGTTAGAACCTACCTCCCTCGCCGGGCA...............................................
...............................**********************************************...............................................
...............................**********************************************........................................A......
...............................**********************************************...........................................T...
,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
                                                                             ...............................................
                                                                             ...............................................
                                                                             ...............................................
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

At first glance there appears to be only a single read spanning this insertion, although the abrupt end of many reads at the insertion breakpoints made me suspicious. Indeed, upon closer inspection, 19 additional reads have primary alignments with soft clipping at the insertion breakpoints.
By default, soft-clipped portions of the reads are not displayed. Is there any way to toggle this behavior?

How would they be displayed? There's nothing in the CIGAR strings to align the soft-clipped bits to the reference

They could be displayed in the same way that igv does it: as a string of unassigned bases at the end of each read.

A program called gridss does this. It can be found here. It basically extracts the clipped bases and repeatedly realigns them to reference with the aligner of your choice with bwa as default. It will do this until there are no further realignments. The simple command-line for this tool after install is:
java -Xmx512M -cp gridss-VERSION-with-dependencies.jar gridss.SoftClipsToSplitReads I=your_input.bam O=your_output.bam REFERENCE_SEQUENCE=your_reference.fa

I hope this helps!

Could you please explain why you downvoted?

I haven't voted up yet since I haven't had a chance to try out gridss. But yeah, that's an unfair downvote.

This answer does not appear to answer the OP's question?

@winni2k would it be more beneficial for me to have posted an answer that says, "No"? OP's question boiled down to, "Is there any way to toggle this behavior?" There is no way to toggle this behavior, so I proposed an alternate program as a solution. Seems reasonable to me.

Just explaining why someone might down vote. I think it's great that you tried to help though.


text:
How to get log2 fold change of RNA-Seq data for time series experiment?
I know if there is one control and one treatment group it is pretty straight forward to interpret log 2 fold change. But, I have time course experiment. I have infected cells with viruses and I harvested cells at different time points: o hr, 6 hour, 12 hr and 24 hour. Here o hr is my non-infected control. What is the best way to show log 2 fold change for each gene in this type of time course experiment. I ultimate aim to find highly expressed genes at different time points.

If you just want the change with time then use time as a continuous variable. The log2FC will then be the change per hour (you can use the default Wald test to extract it).
For time-point specific comparisons you'll instead need to treat the time points as groups and use contrasts as appropriate.


text:
Adding a string to start of one column
I have this dataframe
> head(mut)
  sample chr     pos ref alt
1  t_005   1  928006   G   A
2  t_005   1 1649842   G   T
3  t_005   1 2020408   G   A
4  t_005   1 2031677   T   A
5  t_005   1 2901878   G   A
6  t_005   1 3039338   G   A
>

How I can add "chr" to the chr column? I mean, instead of "1" I want "chr1", etc.

This is basic R question and more suitable to StackOverflow main site, and it is already has been asked many times, see here for example

This should work:
mut$chr <- paste("chr", mut$chr, sep = "")

if you use paste0 instead of paste, you can omit the 'sep=""'


text:
Plot a circos plot to show the consistency between 2 samples
I have called SNP and INDEL in two matched samples by strelka and extract this information from .vcf file and I have these
CHROMOSOME  POS REF ALT SAMPLE
1   928006  G   A   t_005
1   1649842 G   T   t_005
1   2020408 G   A   t_005
1   2031677 T   A   t_005

and
CHROMODOME  POSITION    REF ALT SAMPLE
1   14115   A   T   o_005
1   541052  T   C   o_005
1   1088123 T   G   o_005
1   1232501 A   G   o_005

I need to show if and how extend the mutations between these two samples are consistent by a circos plot but really I don't know how people do that
Something like this

or

FWIW, I've never found a use case where a Circos plot was the best visualisation tool. They look pretty, but they don't really convey much information.

Thank you, I tried a waterfall to show the landscape of mutations in two samples but I failed again. I need a way to show which extend the mutations between these two samples are consistent

Then you could use a confusion matrix, e.g. http://i66.tinypic.com/eanz9c.jpg

What I should put in confusion matrix? For example total number of mutation? But maybe the total number of mutations is the same but the mutations are different  :(

nonono. Every variants gets one position in the matrix: in my example, a variant which is heterozygous in Sniffles and homozygous in the truth set ends up in the square with 516. A perfect match has only variants on the diagonal.

Sorry do you know any blog of tutorial to showing use case?

I've posted my suggestion as an answer to not clutter the comments here.

I agree that maybe a confusion matrix would be better, however there is a nice circos package in R https://github.com/jokergoo/circlize. You can use the spread function in dplyr to change the dataframe into a square matrix.

Sorry I could not find any vignette for that to adapt my data

A circos plot is most likely not the most appropriate solution here. What I would suggest is a confusion matrix, of which you can find an example here:

For every variant in your vcf you'll add a number in this matrix. One sample is the columns, the other is the lines. If your variant is homozygous in both, then you add in that square +1 (the cell with 5845 in the example).
A perfect concordant sample pair will have only variants on the diagonal.
Here is some python code to get such a matrix. It uses cyvcf2 and pandas, and expects as input a vcf file with both samples.
from argparse import ArgumentParser
from cyvcf2 import VCF
import pandas as pd

def main():
    args = get_args()
    confusion_matrix(args.vcf)

def confusion_matrix(vcff):
    """
    First level of the dict is the "first" call, second level is the "second" sample
    0: hom_ref
    1: heterozygous
    2: unknown/nocall
    3: hom_alt
    """
    zygosities = {0: {0: 0, 1: 0, 2: 0, 3: 0},
                  1: {0: 0, 1: 0, 2: 0, 3: 0},
                  2: {0: 0, 1: 0, 2: 0, 3: 0},
                  3: {0: 0, 1: 0, 2: 0, 3: 0},
                  }
    for v in VCF(vcff):
        zygosities[v.gt_types[0]][v.gt_types[1]] += 1
    zygs = [2, 0, 1, 3]
    df = pd.DataFrame(index=zygs, columns=zygs)
    for tr in zygs:
        for te in zygs:
            df.loc[tr, te] = zygosities[tr][te]
    df.columns = ['nocall', 'hom_ref', 'het', 'hom_alt']
    df.index = ['nocall', 'hom_ref', 'het', 'hom_alt']
    print(df)

def get_args():
    parser = ArgumentParser(description="Create confusion matrix of SNV calls")
    parser.add_argument("vcf", help="vcf containing two samples")
    return parser.parse_args()

if __name__ == '__main__':
    main()

Thank you, for example if I have filtered t_005.vcf anf o_005.vcf files in a directory, how this code locate them?

It needs a single vcf with both samples. For that you'll need bcftools merge to combine your vcf files.


text:
ExpressionSet with a custom gene list
After an ExpressionSet is established,
I want it to contain only the genes of my choice ("selected.txt")
What is the best approach??
Thanks

Welcome to tthe site. Could you show us which format is your selected.txt files? And how is your ExpressionSet object ? In this types of question explaining what do you have and what you have tried make it easier for us to help you.


text:
Why doesn't the mouse GRCm38/mm10 refGene genome annotation file contain non-coding transcripts?
I am new to Bioinformatics and I am exploring the refGene.txt files from the UCSC genome annotation database for several species.
My question concerns the Dec. 2011 (GRCm38/mm10) assembly of the mouse genome. I have seen that the Human one (hg39) contains both coding and non-coding transcripts. But the Mouse (mm10) RefGene.txt only contains coding transcripts.
Why is that?
Source of the file: http://hgdownload.soe.ucsc.edu/goldenPath/mm10/database/refGene.txt.gz

Is there a reason you're using refseq rather than the more complete annotations from Ensembl/Gencode?

Yes, it is an exercise in a course. I am to use bash shell and gawk to extract some information from several refGene.txt files.    The question about the absence of non-coding transcripts in mouse is not part of the exercise, but have puzzled me. I have been looking for info about it, but I haven't a clear explanation.

It wasn't the correct file.
I downloaded the file again from UCSC, compared it with the file that I was using, and they were different. This new file contained both information from coding and non-coding transcripts.
To all the kind people that took interest in my problem, I am sorry.
Edited 30/12/2020:
To add a bit of clarity; the link in the question is the correct link from which to download the correct file. The file with which I was working from the beginning was the wrong one, but I don't really know where I got it from the first time.


text:
How to find Wild Type (WT) protein structure on Protein Data Bank (PDB)?
Currently, I am working on Galactose oxidase enzyme for protein docking and I am totally new to it. I have to study how mutation on specific residue will make mutated galactose oxidase better than the wild type. My questions:

There are so many structure available on PDB, how to know which structure is the so called wild type structure? Is the structure with PDB Code: 1GOG belongs to wild type?
I was requested to extract the ligand(substrate) out of the structure and docked it back. But seem all the structures available on PDB do not complexed with ligand. For what I know, the substrate for galactose oxidase could be D-galactose. Could I use the galactose oxidase structure from pubmed and dock it to the enzyme? Btw, what is purpose of extract the ligand out and dock it back?

Welcome to the site. 1) Did you check out the description of the structures? Of which organism are you interested in the wild type of  the protein? 2) Yes you could. Why do you want to extract the ligand from a structure and dock it back? (This seems homework and could be only to you to practise how to use the software)

Yes, I did read the literature and it mentioned that 1GOF is belong to wild type but interestingly, I found another 2 structures (1GOG, 1GOH) that have the identical sequence as wild type (1GOF). That's why I confused which structure to use, btw the differences between these 3 structures are only their ligands. My question is since their sequences are identical, does it means that I could just use either one of them to perform docking?

All three structures are from the same specie but other species could have a different wild type. There is no rule why a single molecule can't be linked to several others but the structure might be different (or change when docked with other ligands). Yes you could use any structure to perform docking but the suggested structures could be different according to the initial state of the structure


text:
Pattern matching and removing from the name
colnames(countdata) = gsub("Aligned.sortedByCoord.out.bam_rem.bam", "", colnames(countdata))
colnames(countdata)
dput(colnames(countdata))
    c("BAF60A_KD_1_S13", "BAF60A_KD_2_S14", "BAF60A_RA_1_S17", "BAF60A_RA_2_S18",
    "BAF60A_VD3_1_S15", "BAF60A_VD3_2_S16", "BAF60C_KD_1_UI_S1",
    "BAF60C_KD_2_UI_S2", "BAF60C_KD_RA_1_S11", "BAF60C_KD_RA_2_S12",
    "BAF60C_KD_VD3_1_S9", "BAF60C_KD_VD3_2_S10", "VC_1_UI_S3", "VC_2_UI_S4",
    "VC_RA1_S5", "VC_RA_2_S6", "VC_VD3_1_S7", "VC_VD3_2_S8")

After removing the above pattern I want to remove this part from each of the colnames _S[anynumber] so that I can end up like this as an example  BAF60A_KD_1 but not sure how to make that pattern and match the same and remove or replace it with empty space. I tried as below which ends up like all the digits are removed.
colnames(countdata) = gsub("[[:digit:]_S]", "", colnames(countdata))

colnames(countdata)

c("BAFAKD", "BAFAKD", "BAFARA", "BAFARA", "BAFAVD", "BAFAVD",
"BAFCKDUI", "BAFCKDUI", "BAFCKDRA", "BAFCKDRA", "BAFCKDVD", "BAFCKDVD",
"VCUI", "VCUI", "VCRA", "VCRA", "VCVD", "VCVD")

Any suggestion would be helpful

I'm voting to close this question as off-topic because the question is purely about regex, not requiring any biological knowledge. You can remove them with names <- gsub("_S.*", "", colnames(countdata)) followed by removing any number or underscore gsub("[:digit:]|_", "", names)

IIRC gsub has a perl=TRUE option for working with advanced regexes. Unfortunately, the [[:digit:]] style regexes are not supported by the base regex engines in most utilities/distros.

@RamRS what do you mean? [:digit:] is a POSIX character class and should work on the vast majority of regex flavors. This isn't a feature of advanced regexes, but a feature of the most basic of them all: POSIX basic regular expressions. Where have you seen it not work? For example, it works in BRE (printf "foo\n123\n" | grep -P '[[:digit:]]'), ERE (printf "foo\n123\n" | grep -E '[[:digit:]]'), PCRE (printf "foo\n123\n" | grep -P '[[:digit:]]').

@llrs I don't think this is off topic. Bioinformatics involves a lot of basic programming knowledge, and I would say that includes regular expressions. You can do a lot of very productive work in bioinformatics without really knowing the first thing about biology. I would say this comes under the same category as general text parsing questions which we've already decided should be on topic.

@terdon I use homebrew-powered GNU-sed on my MacBook Pro, and it won't s/[:digit:]// without the -r flag. That's the reason I added my disclaimer.

@RamRS That should be s/[[:digit:]]//, which works on GNU sed without the -r. The [:digit:] is the special symbol, you need to put it in a character class ([ ], hence [[:digit:]]) for it to work as a character class. s/[:digit:]// shouldn't (and doesn't, on my GNU sed running on Linux) work even with the -r (which just enables ERE). Even busybox sed (which is as basic as it gets) understands [[:digit:]] without -r.

Apologies, I meant 's/[[:digit:]]//'. Weird that busybox sed works fine but my sed doesn't. I'll check again.

Weird. Seems to work on all my computers now. I distinctly recall avoiding POSIX character classes because they worked inconsistently between my computers.

and it really helped me .even though its really a non-biological or bio info question but these minute things helps me a lot

Without regex, split on "_", drop last item, then paste it back again:
sapply(strsplit(x, "_"), function(i) paste(head(i, -1), collapse = "_"))

Or with regex (stealing the code from tools::file_path_sans_ext function which is used to drop file extensions):
sub("([^_]+)\\_[[:alnum:]]+$", "\\1", x)

# [1] "BAF60A_KD_1"     "BAF60A_KD_2"     "BAF60A_RA_1"     "BAF60A_RA_2"     "BAF60A_VD3_1"    "BAF60A_VD3_2"
# [7] "BAF60C_KD_1_UI"  "BAF60C_KD_2_UI"  "BAF60C_KD_RA_1"  "BAF60C_KD_RA_2"  "BAF60C_KD_VD3_1" "BAF60C_KD_VD3_2"
# [13] "VC_1_UI"         "VC_2_UI"         "VC_RA1"          "VC_RA_2"         "VC_VD3_1"        "VC_VD3_2"

multiple ways of doing wow

@krushnachChandra minor note, you will get a quicker and better answers if you post regex/R posts at StackOverflow.

well actually i got banned from posting it on the R forum ,i guess i asked stupid questions.it says from this account no more questions.But im glad for multiple answers coming out

What do you mean "banned"? This is your profile, right? https://stackoverflow.com/users/6342356/krushnach-chandra

yes i do ask question on this forum and other forum but not in the R section may be bad questions

This is what i see when i try to ask question on that forum "You have reached your question limit
Sorry, we are no longer accepting questions from this account. See the Help Center to learn more."

@krushnachChandra Looks like you have a system ban "for asking poor questions", it is temporary, not sure for how long. Relevant meta post regarding ban: https://meta.stackexchange.com/q/86997/228487

okay no issues..i as quite naive to ask certain question ..


text:
How to convert ChiP-seq track in .bam file format into an array at fixed bin resolution?
Is there a out of the box to convert a .bam file to a .wig file at specific resolution/bin size?
I was thinking to use igvtools, but I am not sure how it works:
igvtools count -z 5 -w 25 -e 250 alignments.bam  alignments.cov.tdf  hg19
What is the diffrence between -w,-z,-e?
-z num  Specifies the maximum zoom level to precompute. The default
               value is 7 and is sufficient for most files. To reduce file
               size at the expense of IGV performance this value can be
               reduced.

 -e, --extFactor num

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

-w, --windowSize num

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

Using deepTools to get a bin size of 10 and 20 threads:
bamCoverage -of bedgraph -bs 10 -b input.bam -p 20

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


text:
Clustering of gene co-expression network by igraph R package
I have constructed a gene co-expression network from RNA-seq data. The network has more than 10000 nodes and more than 1 million edges.
The network file as an edge list format of memory around 1gb which was created by calculating Pearson correlation of each gene pairs and gene pairs having correlation greater than 95% we're selected to create the edge list. The file is too large to run in cytoscape due to limitation of computer configuration.
In R, I am able to calculate the topological properties of the network. And to cluster the gene coexperssion network, I tried using different community detection methods and but I am unable to identify which community detection is good for obtaining clusters?.
After clustering how to rank the subclusters and save each subcluster in separate files.

What kind of clustering do you want? Should the clusters be discrete (each node can only belong to a single cluster) or overlapping (each node can belong to multiple clusters)? You say you want to rank the clusters, OK, but according to what value? Cardinality? Connectedness? Something else?

Why not just plot a direct correlation heat map? Very easy indeed via R

For the first part, do you mean that the file is too large to be run on your computer?
For the second, if I understood correctly, you can use igraph or ggnet2 to color code the nodes based on the conditions you have and see if communities are formed based on condition. This would be a simple way to do this kind of thing.
I personally prefer using ggnet2 as I found it more user friendly. What you can do if you have an adjacency matrix, say adj_mat is to create a color vector based on the conditions for each node in the network, where length(color_vec) = ncol(adj_mat) = nrow(adj_mat) - since it is a square matrix. Following that, you could do something like this:
library(ggplot2)
library(network)

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

color_vec <- as.factor(color_vec)

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

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

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

I think that Cytoscape has problems with dealing with such big networks even if OP had the computational resources  needed.

I don't have much experience with Cytoscape.. Also I have never had to create a network this large. Would it not be feasible in R?

Yes, it is possible to analyse it with igraph. But in which format do you currently have the network and how do you create it?


text:
How to use busco installed through conda?
I have installed busco on my Mac computer using conda install --channel bioconda busco=3.0.2 into an environment called busco3_0_2. However, when I type busco, the command is not found.
Anyone knows what the name of the binary is? Usually when I install something with conda, the name of the command is the same as the name of the package. In this case, it does not seem to be so.
Note that the package installation finished successfully.

From the README, it looks like the command-line entry point is not busco but run_BUSCO.py. This seems to agree with the bioconda recipe, where run_BUSCO.py and run_busco are listed among the commands.

Cheers. I found run_BUSCO.py and run_busco in ~/anaconda/envs/busco3_0_2/bin.

The new version of busco (v4) has the user entry point called busco (it's still a python script). BUSCO v4 is available at conda, note it's a very fresh update.


text:
What does this genome file from NCBI means?
I obtained a genome from NCBI. For instance, Haloquadratum walsbyi strain AM180088.1. But the FASTA file contain two parts. The first part marked as
>AM180088.1 Haloquadratum walsbyi DSM 16790 complete genome
The second part marked as
>AM180089.1 Haloquadratum walsbyi DSM 16790 plasmid PL47 complete genome
What does that mean? In particular, the second one "plasmid"? If it is called complete genome should it all in one? Why divide into parts?

Archea can have plasmids last time I checked...

what does that mean? So, plasmid have different gene from other part of archaea?

Yes, bacteria and sometimes both archea and eukaryotes can have plasmids (google the term if you don't know it) in addition to their normal genome.

I know they have plasmid. But will the plasmid have different gene from the bacteria/archaea itself?

Typically, yes.

I'm voting to close this question as off-topic and not worth migrating to the biology stackexchange.


text:
Viral genome assembly using broad viral ngs pipeline?
I am trying to assemble RNA virus genome using Broad Viral NGS pipeline BROAD VIRAL NGS PIPELINE. I am two questions :
1) As this pipeline requires unaligned bam format as input, how do I convert fastq to bam ?
2) Some of the tools like   assemble_trinity require Trimmomatic clipDb. How can I get this clipDb ?

Since you are already using the Broad Tools sets you can use Picard FastqToSam to make the conversion
As far a clipDb I am unfamiliar with that and a quick google search and look at the trimmomatic manual were unhelpful


text:
taxon exclude list for searching local blast database using blastn
I am looking for a solution to exlude certain entries when searching a local blast nt database (with blastn), specifically the sequences from uncultured / environmental samples, ideally using their taxon ids.
This should ideally replicate the results from the BLAST web interface, when checking the option "Exclude Uncultured/environmental sample sequences", see screenshot below.
Three approaches that I am aware of so far are:
(1) grep taxon IDs for any taxa with "uncultured" or "environmental" in their scientific name from names.dmp and fetch the GIs belonging to these taxa using the E-utils. Then use the -negative_gilist parameter for blastn. However, GIs are discontinued for future records, so this method should be replaced by Acc.Ver identifiers.
(2) locally format the blast database after prepocessing the nt file by myself and exclude records belonging to the taxa list as obtained in (1).
(3) Postprocess the blast output and remove the relevant lines, but this approach might remove a lot of alignments, whereas (1) and (2) would always generate the desired number of output alignments.
So, I am looking for other suggestions, preferably without preprocessing the nt file.
Ideally, one would like to have a -negative_taxalist option for blastn, which would resemble what is possible on the NCBI BLAST web interface, where one can exclude taxa, see screenshot below.

Thanks!

Hurray, the solution seems already to be implemented in the new BLAST database format (v5), which is usable with blast+ versions >= 2.8.0.
It is possible to specify an exclusion as well as inclusion list of taxon ids, which, however, have to be on species level.
Species level taxon ids can be obtained from higher level taxon ids using the script get_species_taxids.sh, e.g.: get_species_taxids.sh -n Enterobacterales
See docs here: https://ftp.ncbi.nlm.nih.gov/blast/db/v5/blastdbv5.pdf
Edit:
For creating taxon ID lists locally from the NCBI taxonomy files, the program subtree can be used. See example therein for compiling a list of taxon IDs from environmental or unclassified taxa.

Sorry to bother you, but how did you manage to work with get_species_taxids.sh? Because I installed edirect on linux, but won't work with the script.


text:
Merging/concatenating index files
I was wondering if any one has had success merging pre-generated reference index files for BWA, bowtie2 from multiple references?
For a concrete example I have hg19 and a bunch of viral and bacterial genomes that all have been index separately but I want to do a single alignment to all of them.  Is there a way I could say combine hg19, virus_A, virus_D and bacteria_X into a single reference without having to regenerate the indexes?  Since the problem could get combinatorial quickly and indexing anything with hg19 is time consuming I was wondering if there was a way around it.

this will likely depend on the tool no? i doubt there is a generic why to merge indexes? for BLAST, there is blastdb_aliastool

Yeah, I would split this into separate questions for each tool. I don't actually know how the blast index works, but I would be very surprised if the same approach will work for all three tools. It might be the same for bowtie and bwa, but blast will likely be different.

Agreed I wasn't looking for a generic tools but the question does read that way.  I was looking for a way for each of those tools.  Thanks @Chris_Rands ill check out blastdb_aliastool

What would you be aligning? Short reads? Why would you want to align the same data against three so different genomes? Surely the time you would save by not re-indexing would be lost again while wasting time aligning to the wrong genome, right?

@terdon I imagine they have metagenomic data or similar where the sample composition is not known a priori (reads could be 90% human, 8% E. coli etc.)

@Chris_Rands duh, indeed, should have thought of that, thanks.

@terdon Correct we don't know the make up of the sample and want to have the reads complete for alignments

That makes perfect sense for your BLAST question but I would be wary of using it for short reads. Even between bacteria and humans there are bound to be regions of homology. I would be very curious to know how well this works if you do manage to do it! Maybe by assigning the viral and bacterial genomes as extra "chromosomes" in a single multi-fasta file along with hg19 and indexing that?

@terdon The multi-fasta is the base case I was hoping not to have to reindex the different combinations but will either do that or create one giant reference

Could some of these help? 1, 2, 3?


text:
Merging / Concatentating BLAST Databases
I was wondering if any one has had success merging BLAST databases from multiple references?
For a concrete example I have hg19 and a bunch of viral and bacterial genomes that all have been index separately but I want to do a single alignment to all of them.  Is there a way I could say combine hg19, virus_A, virus_D and bacteria_X into a single database without having to regenerate?  Since the problem could get combinatorial quickly and making a db with hg19 is time consuming I was wondering if there was a way around it.
This is similar to another ticket but asking for BLAST instead of alignment indexes.

There's an app for that, apparently.

Often we need to search multiple databases together or wish to search a specific subset of sequences within an existing database. At the BLAST search level, we can provide multiple database names to the “-db” parameter, or to provide a GI file specifying the desired subset to the “-gilist” parameter. However for these types of searches, a more convenient way to conduct them is by creating virtual BLAST databases for these. Note: When combining BLAST databases, all the databases must be of the same molecule type. The following examples assume that the two databases as well as the GI file are in the current working directory.
Aggregate existing BLAST databases
To combine the two nematode nucleotide databases, named “nematode_mrna” and “nematode_genomic", we use the following command line:
$ blastdb_aliastool -dblist "nematode_mrna nematode_genomic" -dbtype nucl \
  -out nematode_all -title "Nematode RefSeq mRNA + Genomic

So, you can use the blastdb_aliastool to create merged databases just like you asked for.


text:
How to find the position-specific weight matrix (PSWM) for a list of genes having different sizes?
Wiki has the algorithm for finding position-specific weight matrix (PSWM) for a list of nucleotide sequences having an equal number of nucleotide. How do we find the PSSM for a list of nucleotide sequences having different sizes?

If your sequences are different lengths, this implies a gapped alignment with insertions/deletions. This is not possible with a regular PSSM, so you have to create a separate PSSM for each group of homologous sequences of the same sequence length.
The alternative is to use a different approach that allows gaps like profile HMMs (as implemented in HMMER3 or HH-Suite3).


text:
residual Squared Coefficient of Variation (rCV²) vs Distance to Median (DM)
I am working with single cell RNA-seq data. I obtained the squared coefficient of variation (CV²) as a measure of gene expression variability:

I want a metric to express gene expression noise that eliminates/takes into account the observed relationship between CV² of a gene and its mean expression so I can compare this metric across genes/conditions:

So far I have found two metrics that try to tackle this problem:

Distance to median (DM):

A median-based trend is fitted to the log-transformed CV² against the log-transformed mean. The DM is defined as the residual from the trend for each gene. This statistic is a measure of the relative variability of each gene, after accounting for the empirical mean-variance relationship

residual CV² (rCV²):

Fit between the mean log2 expression and CV² using a gamma generalised linear model as parametrisation then computing the absolute deviation of the CV² values from the fitted CV² values.

Do you have any clue about which metric is better or if there exists any other similar metrics?

There has been a lot of work done on this problem already, particularly these two papers:
https://www.biorxiv.org/content/10.1101/576827v2
https://www.biorxiv.org/content/10.1101/574574v1

You should do both. Distance to median (DM) is the most powerful transformation IMO.
The important feature of the coefficient of variation which is more famously known as the index of dispersion (sd2/mean) is if:

IoD is approximately 1, the distribution is random (Poison),
IoD <1 the distribution is a uniform distribution (very rare in biological systems)
IoD >> 1 the distribution is skewed (clustered), such as a negative binomial distribution. This is common in biological systems.

An example of a uniform distribution are the distance between apple trees in an orchard. In general biological statistics IoD is rarely useful because many distributions are moderately over dispersed. This could be different however for transcriptomics. Thus if IoD results in a moderately skewed distribution you can comment on it, but its not worth taking further.
Distance to median appears to be a "non-parametric" standardisation and that is exactly what you are looking for here. Standardisation is a formal process of accommodating natural variation, e.g. time of day of the experiment, to extract the underlying trend. It is very common (in fact almost essential) in time series.


text:
BUSCO does not find buscos in reference sequence of C. elegans (technical problem)
I am working on a GNU/Linux machine. I have installed BUSCO using conda. I have downloaded the nematode dataset nematoda_odb9.tar.gzhere. I downloaded the reference sequence of C. elegans from Ensembl. I then ran busco using the following command:
run_busco --in Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --cpu 1 --out busco --lineage_path ./nematoda_odb9/ --mode genome --force > busco.o 2> busco.e
tblastn seems to have run properly. run_busco/blast_output/tblastn_busco.tsv contains lots of hits for many of the buscos.
However, I have 982 "Missing BUSCOs".
The file that contains logs and errors related to augustus (run_busco/augustus_output/augustus.log) contains:
Will create parameters for a EUKARYOTIC species!
creating directory ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/ ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_parameters.cfg ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_weightmatrix.txt ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_metapars.cfg ...
The necessary files for training BUSCO_busco_1395218643 have been created.
Now, either run etraining or optimize_parameters.pl with --species=BUSCO_busco_1395218643.
etraining quickly estimates the parameters from a file with training genes.
optimize_augustus.pl alternates running etraining and augustus to find optimal metaparameters.

It does not seem to contain any errors.
I have the following error:
Failed to open sequence file run_busco/augustus_output/extracted_proteins/EOG091H07ZC.faa.1 In fact, the folder run_busco/augustus_output/extracted_proteins/ is empty.
Anyone would know why I have no extracted proteins?

Most likely your Augustus installation is not functioning properly. Did you install Augustus from source or using conda? I believe the conda version of Augustus is known to have some issues.
Something worth trying: Running a complete busco install in a VM. This image has all dependencies installed and properly configured.

It was installed via conda, the log output clearly shows the Anaconda being accessed for busco installation. Personally I'm surprised because Anaconda programs are very reliable, but that might just be the core installation (at least for  Mac OSX). Anyway its good advice to try a direct installation.
Could you please point out what "some issues" might be?

BUSCO worked with this Docker image: https://hub.docker.com/r/comics/busco I'll try again using the conda installation and point towards some issues but for the moment, I am happy I could run it in a Docker container.

@MichaelG. you are right that is quite reliable in general. Augustus however, is quite tedious in it's installation (e.g. the location of the config file and the location of it's species training-sets). For example, the manual of BRAKER(the automated gene annotation program) states: "Problems have been reported from users that tried to run BRAKER with AUGUSTUS releases maintained by third parties, i.e. Bioconda". They therefore recommend to install Augustus as distributed by the original developers

@charlesdarwin That's good to hear! Most likely the person that build that docker container did a correct installation of Augustus for you :)


text:
Is the visual cortex of a newborn baby immediately capable of object detection or is this skill learned over time, and if so, how?
Is the visual cortex of newborn babies right off the bat capable of making sense of raw visual data, for instance, converting the constant stream of raw RGB images perceived by the eyes into a meaningful higher-level representation of objects in motion in a 3D world? Yes? No? If not, then does it mean this skill has to be developed over time by means of some learning mechanism?
If the visual cortex needs time to learn advanced visual skills, how does the actual learning mechanism work? Does the visual cortex have to optimize the connections between neurons, in a way analogous to how artificial neural networks optimize their parameters through backpropagation algorithm (machine learning)? If this is the case, then where does the visual cortex get its error signals from? To make the last question clear, in machine learning the typical approach is to compute the gradient of a loss function which compares the model's prediction with the ground truth, and the model's parameters are updated by moving the parameters in the direction of the gradient. If the visual cortex is learning advanced visual skills by virtue of a similar learning mechanism, then what kind of loss function is the visual cortex optimizing?

This isn't on topic here (I am trying to find a better home for it) but what in the world gave you the idea that animals don't need this sort of training? That's precisely what games like this are for, for example. And that's one of the benefits of play in young mammals of all species.

I'm voting to close this question as off-topic because it does not appear to be about bioinformatics. Maybe biology.SE?

The question relates to developing better deep learning (neural network [NN]) machinery based on human/animal perception for image recognition (convolution is a classic NN method). This idea has allegedly succeeded before in NN. However, the empirical principles of why a DL NN approach works are often poorly understood, the OP will be performing vast numbers of differential equations, in contrast biological approaches are better understood, so I'm not sure this question is going anywhere in any case.

Just to finally add that the modern rule of thumb in image recognition is that the deeper the NN the better the result. So just increase the depth of the NN by ten-fold.

It is a good question, interesting too. Now I see that motivation is very computational, but all the questions are very biological and not sure if there is anybody here with the expertise to answer half of them.

Also, it's super broad. It would be a decent thesis topic.

While the question as is is now more about bioinformatics I think it is too broad to be answerable in this format (we can't answer a thesis in a single answer). Perhaps you have a particular question you are having problems with (this seems like too general to be an actual problem you have)


text:
How to allow "half windows" in biopython?
I am using Biopython to examine protein sequences including the Kyte and Doolittle hydrophobicity. However, for short sequences, this returns an empty value; the value is not allowed if the window does not have 20 values.
from Bio.SeqUtils.ProtParam import ProteinAnalysis

window_length=20
edge=1
kyte = {'A': 1.8, 'R': -4.5, 'N': -3.5, 'D': -3.5, 'C': 2.5, 'Q': -3.5, 'E': -3.5, 'G': -0.4, 'H': -3.2, 'I': 4.5,
        'L': 3.8, 'K': -3.9, 'M': 1.9, 'F': 2.8, 'P': -1.6, 'S': -0.8, 'T': -0.7, 'W': -0.9, 'Y': -1.3, 'V': 4.2}

full_sequence_analysis = ProteinAnalysis("CYGTGAVYLLGLGIGGFS")
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)

Is there are a built in Biopython option in this for "allow half windows/short windows" for incomplete windows? There is nothing I can see to suggest so in the docs.

@Chris_Rands This would fix the error. However, it would change the absolute values. Ideally, I would like to keep the window length at ±10 residues (up to 20 residues) if possible.

Have you tried odd values for window_length? https://web.expasy.org/protscale/ allows 3 5 7 ... 21

@nicolallias The issue would still remain. Even if I had 3 residue window length when I sample the final residue in the sequence, there are only 2 residues available for the window (the final residue, and the preceding residue), so no value is generated because 3 residues are needed for the windowed value.

You can choose between a window_size of 20 for sequences longer than 20 AAs, and the length of the sequence for those shorter::
seq = "CYGTGAVYLLGLGIGGFS"
window_length = min(len(seq), 20)

full_sequence_analysis = ProteinAnalysis(seq)
kyte_window = full_sequence_analysis.protein_scale(kyte, window_length, edge)


text:
Executing PyMOL from a Shell script
I'm trying to execute pymol from a shell script (and it is not working). I'm not executing the script on  PyMOL but even if I do this, it doesn't work neither. However, if instead of a script the commands are manually written in PyMOL, everything goes well.
Here is the code:
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
' > scriptt.pml
pymol scriptt.pml

One of the problems, you will see, is that PyMOL does not remove the water molecules. Another issue is that the program does not exit. Do you have any idea what's going on there?

What program does not exist? What errors are you getting? Is pymol in your PATH? I have no experience with pymol, are you sure remove resn HOH is the right command to remove water?

@terdon the post reads 'exit', not 'exist'. And the same commands typed in the program interface, works fine...

I'm thought Pymol does hydrate and and remove water. I dunno about the command line, but there's options in the GUI

Your shell script is fine, and works as expected. The reason it hangs is that you aren't telling PyMol to exit after running your script, so it's still running and waiting for the next command. Just add a quit at the end and it should work:
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
select prot, 1OYT and not ligand,
save lig.pdb, ligand
save prot.pdb, prot
quit
' > scriptt.pml
pymol scriptt.pml
obabel prot.pdb -xr -O prot.pdbqt
obabel lig.pdb -O lig.pdbqt
echo 'babel conversion is finished'


text:
The best tool to count the number of insertion sequences in a set of bacterial genome
I have a set of 414 GenBank genomes and I need to count a number of insertion sequences (ISs, transposases) in each of them. Is there any traditional way to perform such an analysis?
I came up with three possible ways:
1) Count a number of "transposase" as a product of CDS for all CDS in annotation file of each genome. -- it depends on annotation
2) Use a tool named OASIS. -- it depends on annotation
3) Use a web-server named ISsaga. -- it can analyse only 10 replicons at once
4) Perform blast(n or p) against ISfinder database and count a number of hits with low e-value. -- it seems not reliable
Any help is appreciated!


text:
Access base aligned to particular reference position
The short version: If I have a SAM record, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
The long version: I'm using pysam to write some genotyping code. I have a BAM file with reads aligned to some amplicons. I'd now like to genotype a couple of positions of interest in each read while maintaining the phase information. For example, at amplicon 1, if read 1 has A for SNP1 and T at SNP2, then I want to report an A,T genotype for read 1 at amplicon 1.
The typical way to genotype SNPs of interest is to compute a pileup and look at all of the aligned bases position-by-position. The pysam API includes the following example of how to work with pileups.
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    print ("\ncoverage at base %s = %s" %
           (pileupcolumn.pos, pileupcolumn.n))
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            print ('\tbase in read %s = %s' %
                  (pileupread.alignment.query_name,
                   pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

However, I lose the phasing information if I analyze the data position-by-position. Is there any way to flip this example so that I'm iterating over reads in the outer loop, and then over positions in the inner loop? In other words, if I'm iterating over reads, is there any simple way to retrieve the base aligned to a particular reference position without computing a pileup?
As a potential solution, I've though about storing genotype information in a nested dictionary indexed first by read ID and then by position, but depending on coverage that could be incredibly memory intensive.

You can use get_aligned_pairs(); you'll get a list of tuples of (position_in_query_sequence, position_in_reference_sequence)

That's promising! It's been a hot second since I've worked on this code, but it looks like it could indeed allow me to keep phase information on a read-by-read basis. I'll try it out when I can.

You can use the get_reference_positions() function in pysam to get a vector of the positions. Searching that for your positions will allow you to output phased calls.
for read in bam.fetch(some_position):
    positions = read.get_reference_positions():
    if some_other_position in positions:
        idx = positions.index(some_position)
        idx2 = positions.index(some_other_positions)
        ... and so on...

That's not exactly correct syntax, but you get the idea.

Thanks, this is helpful. When there are no insertions (I blocks) in the CIGAR, I can use idx and idx2 to index directly in to record.query_sequence and get the base at the correct position. However, reads with insertions need additional operations to line up the get_reference_positions() indices with the query_sequence indices. So unless I'm missing something, it looks like there's a way to do this without computing a pileup...just not a simple way. :-)

I think you summarized it nicely, possible but not simple.

I had the same problem as you, and so I wrote my own code and sharing it here. No guarantee that it's bug-free ;-)
This code takes advantage of the cigar string. It returns None if the base has been deleted or clipped, or if the bam file did not store the sequence.
from typing import Optional
import pysam

def find_base_in_alignment(alignment: pysam.AlignedSegment,
                           pos: int,
                           bam_stores_revcomp: bool = False) \
        -> Optional[str]:
    idx_q = 0
    idx_r = pos - alignment.reference_start
    if bam_stores_revcomp:
        seq = alignment.query_sequence
    else:
        seq = alignment.get_forward_sequence()

    if seq is None:
        return None

    for op, l in alignment.cigartuples:
        ref_consumed = op in {0, 2, 3, 7, 8}
        query_consumed = op in {0, 1, 4, 7, 8}

        if ref_consumed:
            idx_r -= l
        if query_consumed:
            idx_q += l

        if idx_r < 0:
            if query_consumed:
                # base is in query between idx_q-l , idx_q
                base = seq[idx_q + idx_r - 1]
                return base
            else:
                # position has been deleted
                return None

Let's say you want to genotype all reads by looking at the base at chromosome chrom and position pos:
bamfile = pysam.AlignmentFile(filename, "r" )
for alignment in bamfile.fetch(chrom, pos-1, pos+1):
    base = find_base_in_alignment(alignment, pos)
    # ... do something with it

Note that the parameter bam_stores_revcomp is important. Some aligners will store the original query sequence in the bam, and other aligners (like minimap2) will store the forward sequence (i.e. the reverse complement whenever a read was aligned to the reverse strand). If you use the former, set bam_stores_revcomp to False, otherwise set it to True

Hi, thank you for your code! I think it should be if idx_r <= 0: , then your code would report the exact same data as IGV, excluding insertions. Thank you!

this should be the accepted answer.

This is not [yet] a complete answer, but hopefully it helps get you on the right path.
I wrote some code to retrieve the subsequence from reads in a BAM file that are actually mapped to the reference, which required calculating mapped base locations. Here's a subsection of the code that does this (slightly paraphrased):
  my $refPos = $F[3];
  my $cigar = $F[5];
  $cigar =~ s/[0-9]S$//;
  my $seq = $F[9];
  my $qual = $F[10];
  my $startTrim = 0;
  my $matchLen = 0;
  while($cigar =~ s/^([0-9]+)([MIDNSHP=X])//){
    my $subLen = $1;
    my $op = $2;
    if($op eq "S"){
        $seq = substr($seq, $subLen);
        $qual = substr($qual, $subLen);
      }
    }
    if($op =~ /[M=XI]/){
      $matchLen += $subLen;
    }
  }
  $seq = substr($seq, 0, $matchLen + $addSeq);
  $qual = substr($qual, 0, $matchLen + $addSeq);

So, some things to take note of:

The entire CIGAR string (excluding terminal soft clips) should be parsed up to the required base location
Any soft clips at the start of the sequence need to be accounted for
Only update the read location counter if there's a match, mismatch, or insertion


text:
How to filter blast results from blast run against virulence factor database
I  downloaded virulence factor database from here , to predict virulence genes from few genomes that my lab  sequenced. I  performed local blast on the database with default setting and found around 4000 genes. I want to filter this result. I am not sure how to proceed. There are evalue, pident (percent identity), qcovs (query coverage) etc that I can tweak and I did tweak. For example, if I set pident above 40% the qcovs automatically goes close to 100%. That gave me indication above 40% is pretty significant and biologically relevant. But I am not sure what percent I should consider for treating the query as having the same function as the subject -- in this case that function is virulence.
Need some guidance.
Edit1: I ran blastp
Edit2: I am blasting protein (cds predicted from genome by annotation tool prokka). My question here is what is the most reasonable way to filter putative virulence factor using blast.

Please edit your question and tell us what your queries are. What are you blasting? Proteins? DNA? cDNA? Are your queries known genes? Are they the sequences from the virulence database? Why are you using blast for this? What is the question you are trying to answer?

I am blasting protein (cds predicted from genome by annotation tool prokka).

And how close are the genomes you are blasting? You can't have an absolute value to filter by, it will always depend on the phylogenetic distance between query and subject. How many proteins are you comparing? And why blastp if you are blasting against genomes? Do you mean tblastx or blastx? Also note that virulence is not a function. What you want to do, presumably, is identify the homologous proteins. Blast isn't necessarily the best tool for this. Have a look at my answers here and here.

I am using local blast and blasting against virulence factor database found here http://www.mgc.ac.cn/VFs/main.htm.

"Do you mean tblastx or blastx"
No blastp. I have extracted translated cds from annotated gbk files.

The answer is it depends. 40% pident with ~100% qcovs (assuming length of protein is not very small) tells you that the virulence factor (VF) and your prokka predicted protein are likely homologs, they belong to the same (or closely related) protein family. For some research studies, you might annotate your prokka proteins as 'putative' VFs.
On the other hand, in a clinical context, even 90% pident might not be enough evidence to say these predictions are real VFs. A single amino acid change alters the protein and therefore potentially the function. You could investigate bioinformtiaclly how the changes impact the protein structure and domains etc. Ideally validate it experimentally, e.g. for a new putative E. coli toxin, clone the protein and express it in a commensal E. coli strain and test the toxicity of the new strain.
In practise, think about what question you're asking and who you want to convince. Try numerous thresholds and see how robust your results are to threshold changes (4,000 VFs seems like a lot for a 'few' genomes). You could also try the VirulenceFinder database. When annotating antibiotic resistance genes, if you use 80% pident and qcovs thresholds with a minimum protein length then most researchers will accept your predicted annotations as 'conservative', maybe the VF community has similar arbitrary but pragmatic standards.


text:
Viral genome finishing
I have assembled poxvirus genome using Ray. The assembly is good. Out of several thousand contigs I got, I was able to get one scaffold using Contiguator tool, which is about 90% of my genome. I have well annotated reference sequence. If there any good tool for filling gaps and checking assembly quality?
Things I've tried so far:

IMAGE from PAGIT. It asks for ace file which I guess is obtained from Newbler. So, I am not able to use this.
gapfill_py script from broad: There are some gaps still missing.

I came to know that Consed is good. But, I have windows which is not compatible for Consed.
I checked my BAM file with Tablet. There are reads that cover almost entire gene. Are there any good tools or manual ways to finish my assembly so that I can get the full genome?

Is gap filling for viruses in principle different than gap filling in pro/eukaryotes? I could make some suggestions but they are based on my experience in insect genomics.

I don't think it's too different. It's pretty straightforward actually. Virus that I am sequencing does not contain any introns and overlapping ORFs. It has distinct conding regions flanked by promoter sequence in the intergenic region. And I have well annotated reference sequence too. Your suggestion would be really helpful.

Assuming that gap filling is a domain-independent and works the same for viruses as for bacteria or eukarytes.
There are two real options I am aware of: Sealer a successor of GapFiller from ABySS package and GapCloser from SOAPdenovo2 package. According to Sealer paper it's better, but I am not sure if any indepented benchmarking have confirmed their conclusion. On my insect data GapCloser worked better than GapFiller (back then I have not tried to compare it to Sealer).

Can you please tell what input files do Sealer and GapFiller take ?

I believe that in both cases it will be a gapped assembly (contigs glued by NNNNNs into scaffolds) and the sequencing reads that are mapped on the assembly to gradually reconstruct what is in the gaps.


text:
How I debug this code?
I am trying to run CMSclassifier::classifyCMSfunction on my data but I am getting this error
library(CMSclassifier)

> Rfcms <- CMSclassifier::classifyCMS(my_data,method="RF")[[3]]
Error in match.names(clabs, names(xi)) :
  names do not match previous names

This code classifies gene expression data
Code
But on example data code works
I have attached my data and example data here
Could somebody please help me in solving this error?
My data
Example data

Interesting question. Its a random forest (ML) algorithm. In truth I'd look to directly code this algorithm via Sci-Kit learn rather than use a forumulation in an R library. Just my thoughts.

How did you do to debug this error? Did you try with the example data and you got the same result? Did you try to find what is causing the trouble on the example data? Did you searched how to debug this ? Have you used traceback() to locate the source of the error?

I have done all but very complicated

The problem was I should use ENTREZ rather than gene symbols

Geee, I don't know ... your code and approach looks very solid to me. I thought maybe the classification was using a different name to the training/test set.

I see what you mean but I'm not sure that is the source of the error. I kinda suspect that the example data set you provided works fine. What I think is happening is the labels between the classifier, test and training data are different at some point and so classification ceases.
Your code (albeit I rarely use R) and approach and data look fine. In fact very professional, the data appears scaled (not shown). It looks good.
Perhaps what you refer to is the row name consistency in labelling against the example data set. For example row name, ARHGEF19, this is a gene ID rather than a NCBI code, viz.
efetch -db nucleotide -format gpc -id "LS482351" | xtract -insd gene gene

Output
LS482351.1  ARHGEF19

Elsewhere you use an AL590644.1 as the label, I assume this is the genomic region for PADI4
efetch -db nucleotide -format gpc -id "AL590644" | xtract -insd source organism

Output
AL590644.14 Homo sapiens

(There appears to be version 14 for this subgenomic locus).
So the labelling is a bit mixed and you can feasible convert between them using efetch and scripting a hash/dictionary. However, I don't think its the problem here because its a random forest classification and a label is a label.
It is possible the package you are using is buggy and some of the labels are not getting parsed properly e.g. with a XXX-1 tag or even if a XXX.1 tag. The example data set uses "nice whole" integers. I don't know about R but if regex approaches are used for label parsing things can easily go wrong (its not a good approach). Just a guess however and its easy to check.

If this was Python Sci-Kit learn I'd give the full debugging a whirl. I don't want to get caught up in R machine learning, but recognise you are using R as your ML approach.


text:
How can I ask snakemake to produce a dag where each node represents a rule?
I am using snakemake to create workflows. It is very convenient to visualise my DAG using snakemake --dag target{sampleA,sampleB,sampleC}File | dot -Tpdf > dag.pdf. The resulting pdf shows all the rules' dependencies to get to the target files. However, there is one node per rule/sample combination in the pdf image. When I have many samples, this creates an unreadable image with little added value. Is there an option to merge all the nodes from the same rule together so that I have the workflow for just one sample?

I think you are looking for --rulegraph instead of --dag
  --rulegraph           Do not execute anything and print the dependency graph
                        of rules in the dot language. This will be less
                        crowded than above DAG of jobs, but also show less
                        information. Note that each rule is displayed once,
                        hence the displayed graph will be cyclic if a rule
                        appears in several steps of the workflow. Use this if
                        above option leads to a DAG that is too large.