text:
Help with MinION sequencing data species identification
Hi I'm new to bioinformatics and have just completed my first run on the MinION (long read sequencing Oxford Nanopore Technologies). I was hoping someone could direct me towards R packages, workflow, tutorials or guides that will help me identify species that are present in my sample mainly for fungi as my amplicons had ITS1 and ITS2 regions.
I have been watching Youtube tutorials on RStudio for over a week and I feel like I'm no closer than I was when I started. I have both fastq and fast5 files and I tried to run the fastq files using EPI2ME, however, the database is lacking in fungi. I have tried using NCBI's BLAST but I can only run <5kb reads and most of my sequencing reads are >6kb. Also, I've tried using Galaxy to convert my fastq files to fasta but I'm not sure how to run my long reads for species identification using Galaxy or if that is even possible.
Any help would be greatly appreciated. Thanks.

What sequencing kit did you use? You mention amplicons... are you amplifying specific regions, or doing a whole-genome metagenomic analysis?

There are a number of useful metagenomics tools out there, including Centrifuge and Kraken2, as @gringer has mentioned.
Also, FWIW, the WIMP tool from ONT, a Centrifuge based metageonomic classifier, does include fungus sequences derived from RefSeq, though it may not include the species you are looking for.
If you want to use NCBI BLAST, you could split your long reads into pieces that are <5kb (there are a number of tools out there for FASTQ manipulation - my current favorite tool for this is seqkit, where you could do something like this):
seqkit sliding -s 4900 -W 4900 test1.fastq
The above will split your FASTQ reads into chunks of 4.9kb. You can adjust the -s parameter if you instead want overlapping chunks.
From a GUI standpoint, you can use Galaxy for the above by making a workflow with the following steps:

Trim Sequences (enter reasonable numbers here, to trim each of your reads to <5kb)
FASTQ to FASTA converter
NCBI BLAST+

Another possibility is QIIME2, which can be used with the UNITE fungal training set. See recent reference (looking at the ITS1/ITS2 regions), UNITE, and QIIME2.
Long term, you may want to look into using the command line, as it is much more powerful and customizable, but if BLAST gives you reasonable results for your reads, the above process in Galaxy may work for you.

I've found Centrifuge and Kraken2 to be useful for species classification. There are pre-prepared databases that include Fungi, e.g. see here for the Kraken 2 "PF" database, which is the Standard database plus protozoa & fungi.
Can't really provide much more help than that with more information about what you've done.
You'll need a computer with around 64 GB of memory in order to be able to run these programs, and it's probably not the best idea to jump right into command line compilations when starting out... but it looks like you're up for a challenge given that your first MinION run has involved jumping into the deep end with a fungi metagenomic project.


text:
Pick matching entry from snakemake config table
I'm trying to run a command using snakemake and a TSV config file. Here is my sample config file:
sample_id       FQ1_id  FQ2_id
MyCleanName-1       MyFqName-1_R1.fastq.gz        MyFqName-1_R2.fastq.gz
MyCleanName-2       MyFqName-2_R1.fastq.gz        MyFqName-2_R2.fastq.gz

Here is my snakefile:
import pandas as pd
samples = pd.read_table("samples.tsv").set_index("sample_id", drop = False)

rule all:
    input:
        expand("{sample}_{num}.fastq.gz", sample = samples.sample_id, num = [1,2])

rule my_job:
    input:
        fq1 = expand("{fq1_id}", fq1_id = samples.FQ1_id),
        fq2 = expand("{fq2_id}", fq2_id = samples.FQ2_id),
        logdir = expand("{sample}.logs", sample = samples.sample_id)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
    message:
        "Running MyCmd on {wildcards.sample}"
    shell:
        """
        MyCmd --sample-name {wildcards.sample} -i1 {input.fq1} -i2 {input.fq2}
        """

rule create_logdir:
    output:
        directory("{sample}.logs")
    message:
        "Creating log directories"
    shell:
        "mkdir -pv {wildcards.sample}.logs/{{MyCmd,MyOtherCmd}}"

When I run snakemake -np MyCleanName-1.fastq.gz, I get the following:
Building DAG of jobs...
Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3

[Fri Dec 24 16:38:54 2021]
Job 1: Creating log directories

mkdir -pv MyCleanName-1.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 2: Creating log directories

mkdir -pv MyCleanName-2.logs/{MyCmd,MyOtherCmd}

[Fri Dec 24 16:38:54 2021]
Job 0: Running MyCmd on MDA-IBC-3

        MyCmd         --sample-name MyCleanName-1 -i1 MyFqName-1_R1.fastq.gz MyFqName-2_R1.fastq.gz -i2 MyFqName-1_R2.fastq.gz MyFqName-2_R2.fastq.gz

Job counts:
        count   jobs
        2       create_logdir
        1       my_job
        3
This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

How can I restrict this so the run only creates log files for the Clean Name being used and uses only the corresponding FASTQs? I feel like I'm missing something trivial here. Should I not be using the config file variable samples directly in the rule my_job, but somehow wiring up the corresponding sample_id and FQ_ids on the all rule?

expand is a Snakemake help function that tends to create confusion for beginners. The goal is to generate lists of file names from patterns, usually to generate the input of a "gathering" / "summarizing" rule. If you want a single file name, I would recommend to use standard Python (>=3.6) f-strings.

However, in your case, I think that the use of a function as input will likely be part of the solution: https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html?#functions-as-input-files Maybe you can have a look at this answer for some ideas of how to proceed: https://stackoverflow.com/a/64660447/1878788

Thank you, @bli - I'll look into this - functions as input might make the most sense for my case.

Using pointers provided by @bli, I made a few changes and the Snakefile started working as expected.
Instead of using expand for the rule's input, I switched to a function that accepts the wildcard sample as input and returns a dict with the actual inputs.
 # Could be a one-liner but I wrote it in a verbose
 # manner so it's easier to read
 def get_input_files(in_sample_id):
     fq1 = samples.loc[str(in_sample_id), 'FQ1_id']
     fq2 = samples.loc[str(in_sample_id), 'FQ2_id']
     logdir = "{in_sample}.logs".format(in_sample=str(in_sample_id))
     return_dict = { 'fq1' : fq1, 'fq2' : fq2, 'logdir' : logdir }
     return(return_dict)

And the rule my_job is now changed (just the input part):
rule my_job:
    input:
        unpack(get_input_files)
    output:
        "{sample}_1.fastq.gz",
        "{sample}_2.fastq.gz",
...
...

Hope this helps someone looking for a similar solution in the future!

source: Stack Exchange

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

sort: votes


text:
Best methods for clustering heatmap based on regions where the genes are most highly expressed?
I am analysing RNA-Seq data for a set of genes from different human brain regions using data from the Allen Brain Atlas. I have made a heatmap using GraphPad Prism. In this heatmap, I show the average RPKM value for each gene of interest for each brain region (the average is taken as there are multiple samples for each brain region).
I am new to creating heatmaps but I know that you can use different methods to cluster the heatmap so that the data is more easy to interpret. I want to cluster the heatmap so that the brain regions which have the highest level of gene expression overall are grouped together.
I have seen online that R can also be used to create a heatmap, and I have had previous experience with R. Does anyone have any recommendations on any good tutorials (that are not too difficult) on creating heat-maps in R as well as clustering the data?

For visualization of expression values it usually makes sense to standardize the data first. The problem with plain expression values is that they differ a lot between genes and as such highly-expressed genes dominate the color scale. As you see in your heatmap, the first and last rows barely show any difference because the color scale is dominated by a single gene CHL1. Therefore the clustering would be due to expression level and not due to differences between columns. For this we standardize the data. Assuming your expression matrix had RPKMs (or similar values) we first log2-transform it and then standardize:

mat_log2 <- log2(mat+1)

#/ standardize (Z = Z-scored)
mat_Z <- t(scale(t(mat_log2)))

#/ cluster rows:
hclust_row <- hclust(dist(mat_Z), method="ward.D2")

#/ cluster columns (optional)
hclust_col <- hclust(dist(t(mat_Z)), method="ward.D2")

#/ plot heatmap:
library(ComplexHeatmap)

hm <-
  Heatmap(matrix=mat_Z,
          cluster_rows=hclust_row,
          cluster_columns=hclust_col,
          name="relative expression")
draw(hm)

See the ComplexHeatmap documentation for details on the package itself, it is very powerful but also monolithic. This is just an example, any other heatmap package will do fine as well.
Important is to standardize data to put emphasis on differences between columns rather than on absolute expression levels.

Thanks for your reply. Would it be possible to know why the RPKM values are log2-transformed before the Z-score standardisation is done?

@ceno980 log2 is generally done on genomics data to dampen the effects of outliers as the gene expression values have such a wide dynamic range. As the standardization divides internally by the standard deviation it just helps that single values skew the calculation overly much e.g. in case where some samples have 0 and the others have thousands of counts. The log2 scale tries to make things a bit more homoskedastic. I guess that is more important if you do things like PCA on the expression values directly rather than Z-scoring but for consistency I always do it.

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


text:
Tutorials on phasing and imputing low-coverage sequencing data
I am new to low-pass whole genome sequencing and have the basic idea of phasing and imputation. I have .vcf file after calling haplotypecaller tools from GATK. After searching the phasing and imputation tools I've encountered some tools like Beagle, Minimac, Shapeit, Glimpse, Eagle. The problem is that I didn't find enough hands-on material on how to process a .vcf from unphased to phased and then to an imputed one.
Can anyone suggests any hands-on tutorial regarding phasing and imputation?

What is your sample size, organism, and average coverage per sample?

total sequences read 4093666, Human cell, average coverage (0.1x-6x).

ok, here is the tricky bit: are you trying to impute germline or somatic mutations? I have tried beagle with somatic mutations, and it simply doesn't (and shouldn't) work. If you are using germline mutations, what is your reference panel size?

I am trying to impute germline. I'm trying to impute it with 1000kg phase3 reference panel.

1000 genome phase3 reference panel.

This is all relevant for data 0.5x-1x coverage.
Assuming you have genotype likelihoods data, if you want to phase low-coverage data, the most suitable options are GLIMPSE, Beagle4 and QUILT. Of the three, GLIMSPE and QUILT are the most recent and similar in performance. shapeit, eagle and minimac aren't suitable for low-coverage data.
I have found GLIMPSE to be fast and memory efficient, but it requires several different steps; breaking the genome into chunks, imputing, reforming the chunks and then sampling from the haplotypes to phase the data. It has a fairly good tutorial here.
QUILT has fewer stages and faster, but is unlikely to make much of a difference unless you are working on very large datasets. It has a tutorial here.
Which one you use probably doesn't matter a whole lot, but I would probably suggest using QUILT as it has fewer stages and is slightly more recent.

I am interested in using Beagle. But it doesn't have good documentation and tutorial. Can you suggest any?

Only beagle4 is suitable for low coverage data and it’s pretty old and doesn’t scale well, so I’d advise against using it unless you have a good reason to?

Also, I guess Shapit is only for phasing. I found some tutorial on Shapit+impute2 https://mathgen.stats.ox.ac.uk/impute/impute_v2.html#examples. What do you suggest about it?

Glimpse or Quilt which one is much older?

Also, there is a new Beagle5.2 version. Will it be suitable for low pass sequencing?

GLIMSE is slightly older, Beagle5.2 is not suitable for low coverage and neither is the shapeit2 and impute2 tutorial you linked.

Thanks for your comment. But I found a paper on the performance of phasing and imputation tools https://www.karger.com/Article/Fulltext/489758. The data here has pretty low depth(4x). It suggests that beagle 4 or shapeit+impute2 has a good performance. Can you tell me more elaborately why are you telling that shapeit+impute2 is not a good choice for low pass sequencing? Maybe add your suggestion in the answer?

I suppose it depends on what you consider 'low-pass'. I am referring to data around 0.5-1x coverage, so the results from that paper aren't really relevant to that case, and it predates the release of GLIMPSE and QUILT. Please read the papers I linked for more details. I think that’s everything I’ve got to say about this, as comments are not for extending questions and discussions.

That said, if you want to add in more specific detail about the coverages you are interested in, I will amend the question.

The author of shapeit (Olivier) once told me that shapeit2 imputation simply doesn't work below 4x. I have no reason to doubt that assessment.

I investigated a few different ways of running beagle v3 (same algorithm as v4) for this paper https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4522619/ i found that beagle scaled really poorly in run time with ref panel size. This is why we went just with asian 1kg haplotypes. I think beagle could work for you. There's also quite a bit about that testing in my thesis: https://ora.ox.ac.uk/objects/uuid:19ce0c8e-1d65-44ff-b56c-77ee849b2167

That said, the imputation using beagle on 11k whole genome samples (1.7x) took around 300k cpu hours. I would expect quilt to require less compute time and probably get you similar results.


text:
makeblastdb creating multiple files of unexpectedly large sizes
I have a set of 100 amino acid sequences and I want to perform a BLASTP sesrch against the refseq_protein database. Accordingly I had set up the standalone version of BLAST (Version 2.11.0+) and downloaded the refseq_protein database from NCBI using the following code
wget ftp://ftp.ncbi.nlm.nih.gov/refseq/release/complete/*.faa.gz

The database gets downloaded as 3027 zipped files containing FASTA sequences. I unzipped all these files and concatenated them into a single file refseq_protein.faa (which is around 95 GB in size). Now when I run the following Python code
from Bio.Blast.Applications import NcbimakeblastdbCommandline
from Bio.Blast.Applications import NcbiblastpCommandline

cline = NcbimakeblastdbCommandline(dbtype = "prot", input_file = "D:\\refseq_protein.faa", out
= "refseq_protein")

blastp_cline = NcbiblastpCommandline(query = "D:\\DEP_sequences.fasta', db =
"refseq_protein", evalue = 0.01, outfmt = "7 sseqid evalue qcovs pident")

cline()

response = blastp_cline()

the NcbimakeblastdbCommandline function keeps creating multiple .phr, .pin, .psq etc files which take up a lot of space (In a demo run it had created ~30GB of these files and was still running). I'm afraid this will exhaust the entire space available on my internal hard drive. I'm wondering if there is a way to estimate the total size of the files which NcbimakeblastdbCommandline would create. This will help me in deciding whether or not to switch to an external storage to perform the BLASTP search.
I am aware of the fact that pre-formatted refseq_protein database exists but I'm not sure what value is to be passed in the db parameter of the NcbiblastpCommandline function, because it asks for the name of the database against which the BLASTP search is to be performed. In the approach that I chose, I had the liberty to set the name of the database.
Any suggestions on how to solve this issue would be appreciated.

Welcome to the Bioinformatics site! Hopefully you find the below helpful.
Regarding the creation of a RefSeq protein database locally...as you mentioned, there are pre-made RefSeq BLAST databases available, and looking at the website, the current DB (as of 12/27/2021) is roughly 62.9GB compressed - and is probably ~2x larger decompressed, so I’d guess that the database you are constructing manually will be at least 125GB in size - and likely more.
The above is regarding the technical issue you are currently facing - but I'd back up and look closer at what you are trying to do - align the 100 amino acid sequences to RefSeq. Here are some alternatives that may be -much- easier.

BLAST your sequences on the NCBI website. This may not be as satisfying as building out your local tools and databases, but may get the job done more efficiently. This of course depends on whether this is a one-off thing, or something you plan to do frequently, in which case a local database may make sense.
Download the premade RefSeq BLAST database. You mentioned you were having issues with how to refer to this database (which may be a useful question in and of itself), but generally an alias file exists within the premade database (entitled refseq_protein.pal), and calling this file will allow you to BLAST all the separate chunks. So, the command line might be something like this:

blastp -db refseq_protein -query amino_acid_100.fa

There are a ton more options you may want to use (try blastp -help) , but the above should work.

Use the -remote option in blastp to run from the command line against the NCBI server. The below line should perform the BLAST alignment against the NCBI refseq protein database.

blastp -remote -db refseq_protein -query amino_acid_100.fa


text:
Pros and cons between genomes and transcriptomes?
What are the major differences between the use cases of a shot-gun genome and an mRNA capture transcriptome? Especially when it comes to downstream analyses such as looking for orthology, and selection analysis.
From what I can visualize in my head, I can come up with the following points:

The genome has the structural data, information on chromosomes, and haplotypes while the transcriptome lacks this information.

The genome records the SNPs more directly than transcriptomes. I might be thinking of this wrong but the genomes would have the actual bases and we can then compare the bases between the genomes but for transcriptomes differences in expression would show up as presence/absence variations that have no way of being traced back?

The transcriptomes should be easier to work with when it comes to computational requirements since they tend to be 25 to 50 times smaller than genomes.

I could be completely wrong about the points I listed and could be missing something major. What do you guys think?

It really depends on your use case. If this is human data, where we already have a fairly polished reference genome, then genome sequencing adds value for the genetic variation, e.g. rare disease genetics, GWAS or cancer genomics. Transcriptome sequencing allow exploration of tissue (or cell type) specific patterns and RNA biomarkers etc. Calling variants from RNA can be interesting in the context of splicing and RNA editing.
For a non-model organism where there isn't currently a reference genome, then the genome creates a reference for future studies to map genomic/transcriptomic reads to, so it's usual to do the genome first.
Since you mention ortholgy and selection studies, well traditional comparative genomics studies have done these analyses from the genome only, but the transcriptome can add value depending on your questions.
Finally, I wouldn't consider computational requirements as you stated: it depends on sequencing depth/coverage and other factors, and is likely not to drive the choice between the two assays

Thanks for responding. So, when I mention comparative genomics what I am trying to say is that getting a good genome for a non-model organism is not easy but it seems RNA data is more abundant. I am trying to see if I can stockpile information on plant evolution from RNA data only since genomes are so rare. As for the computational part, when you feed multiple whole genomes to off-the-shelf software they usually just end up with segfaults (since whole genomes can easily be significantly larger than transcriptomes). I am sure that's not an issue in targeted research based on individual samples.


text:
How Does Cox Proportional Hazard Model Describe Interaction between Covariates
I have read STDHA and the documentation for the coxph() function in R. I understand that the cox PH model explores how different covariates simultaneously impact survival. Based on the Lung data example in STDHA, a covariate can seem significantly influential in univariate analysis, but not in multivariate analysis. Following the lung data example, "age" is no longer significant in the multivariate cox model.  Quite literally, this means that age does not impact survival when you account for sex and ph.ecog. However, how does the cox PH model evaluate and decide this? Under the hood, does the cox ph model stratify groups by sex, then evaluate if age is still impactful on survival (and then repeat for ph.ecog)? My intuition says no, because I don't think you can do that with continuous variables (so, in turn, how would the model stratify by age and then test if sex still impacts survival)? I don't need a deeply mathematical explanation, as I'm just a coder!
Ultimately, I want to understand if I can run coxph() on the female and male group separately, or if these results would be redundant since we have already added "sex" as a covariate when analyzing all patients.
This is a conceptual question, but please let me know if I can provide any code. If it means anything, I am using surv & survminer in R.

The simple answer is that Cox PH is a flavor of multiple regression. Multiple regression does indeed try to subtract off the estimated effect of each covariate and look for associations of all other independent variables. Whether this works well in practice is anyone's guess. For a somewhat philosophical high-level discussion see here.
One notable issue is that "significant" is not necessarily a great way to think about statistics, according to professional statistical societies. So, like, if a p-value goes up or down a little bit that's kind of expected when you adjust for other variables, it's not really notable. For one thing, it's closely related to questions of power, driven largely by your sample size.
In simple terms, I would urge you to put all of your data in the model, find some set of independent variables and covariates that you think makes sense before you start analyzing the data, and then build your model accordingly. I don't think that anyone objects to e.g. plotting raw data or model residuals separately in males and females to demonstrate some point, but most associations that are worth caring about will be obvious in a combined dataset.
You can try to explicitly test for an interaction effect in the model. This can be interesting if there's some reason to expect it to be the case, or if there is very compelling evidence that there exists some interaction, but it has the negative of adding to model complexity and further reducing power.

"Under the hood" a Cox model doesn't stratify, consistent with your intuition. It iteratively solves a partial-likelihood vector score function as a function of regression-coefficient values.
It tries to minimize the sum, over all event times, of the difference between each covariate value for the case having the event and a risk-weighted average of the values for all at risk at that event time. The "risk-weighting" is a function of all current regression-coefficient estimates within an iteration. Iteration continues until regression-coefficient estimates are effectively constant. This nicely handles predictors of all types you might use in a multiple-regression model, including interaction terms, while also allowing for individuals that drop out of the study before experiencing the event.
There thus is nothing to be gained by doing separate models for different groups, as Maximilian Press also says in another answer. Stratified models are often used if the proportional hazards assumption doesn't hold, but those model all the data together and often are structured to provide the same regression-coefficient estimates for all strata while allowing for different baseline survival curves among strata.
I also recommend including as many predictors as possible that you might expect to be associated with outcome in your model, but with one caveat. You risk overfitting if you have fewer than 10-20 events per unpenalized predictor in your model (including levels beyond the first of categorical predictors and any interaction or nonlinear terms). So you might have to make some careful choices or use a penalized approach like ridge regression.
For further information, there are nearly 900 answered questions about Cox models on the Cross Validated site.


text:
Nextflow HISAT2 command exit status 255 with prebuild index
I have a problem with my nextflow pipeline.
My workflow looks like this:

Uncompress the hisat2 genotpye index grch38_genome.tar.gz from https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz

Align .fastq.gz files from /https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806626 and https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806627

However, in process HISAT2 I get command exit status 255:

Command error: "grch38" does not exist

The command executed is:

hisat2 --very-sensitive --no-spliced-alignment -x grch38 -U
SRR1806626.fastq.gz > SRR1806626.fastq.sam

Can someone please explain to me what I'm doing wrong?
My script.nf file for reproducing the error looks like this:
#!/usr/bin/env nextflow

nextflow.enable.dsl=2

params.reads = "$baseDir/data/*{1,2}*.fastq.gz" //https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806626, https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR1806627
params.hisat2_index = "$baseDir/data/*_genome.tar.gz" // https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz

log.info """\
         reads:        ${params.reads}
         hisat2_index: ${params.hisat2_index}
         """
         .stripIndent()

process UNCOMPRESS_GENOTYPE_INDEX {

  tag "$hisat2_index_ch.baseName"

  input:
  path(hisat2_index_ch)

  output:
  path('*')

  script:
  """
  tar xvzf $hisat2_index_ch
  """

}

process HISAT2 {

  tag "$reads.baseName"

  input:
  path(hisat2_indices)
  path(reads)

  output:
  path('*.sam')

  script:
  """
  hisat2 --very-sensitive \
         --no-spliced-alignment \
         -x $hisat2_indices \
         -U $reads \
         > ${reads.baseName}.sam
  """

}

workflow {

  Channel
    .fromPath( params.reads )
    .ifEmpty { exit 1, "reads was empty - no input files supplied" }
    .set { reads_ch }

  Channel.fromPath( params.hisat2_index )
    .ifEmpty { exit 1, "hisat2_index was empty - no input file supplied" }
    .set { hisat2_index_ch }

  UNCOMPRESS_GENOTYPE_INDEX(hisat2_index_ch)

  HISAT2(UNCOMPRESS_GENOTYPE_INDEX.out, reads_ch)

  //Command error:
  //(ERR): "grch38" does not exist
  //Exiting now ...

}

Kind regards,
Thomas

Hi @ATpoint, thank you for your answer. You are correct. I changed  -x $hisat2_indices to -x "${hisat2_indices}/genome" and then it worked.

The problem is probably that the tar process creates a folder with the genome files but the input for -x is not the folder with the index files but the basename of the index files itself.


text:
What is the most appropriate way to find the most recent common ancestor between two distantly related species
I want to specifically find the common ancestor between a lobster and a humans. I suspect it was an aquatic worm of some description. But I want to know about the nervous system of this common ancestor.

can you give us some sense of what you've tried? what is the phylogeny you're using, if any? what's the level of resolution you want? For example, if you look e.g. at this little tree (https://www.researchgate.net/profile/Aiasha-Siddika-Arshi/publication/314095464/figure/fig1/AS:466667777073153@1488273676923/Phylogenetic-tree-of-the-Animal-kingdom-2.png), you can find an informatively labeled node, suggesting that the answer might be "the coelomate common ancestor". Another way would be to find an animal phylogeny, read it into R with read.tree, and then use mrca(). IDK which way you want.

Sorry I’m a total novice at this I’ll digest your response and get back to you later.

I guess, it would help to know exactly what you're trying to do. If you're looking for a single piece of factual information, there is a knowable right answer here that you can just look up (which I believe is "the coelomate common ancestor"). However, the fashion that you phrase this and the forum that you post it suggest that you want a more general workflow.

for example, see this web app which computes common ancestors: http://www.timetree.org/. To state it a different way, your specific lobster/human question is finding the common ancestor of deuterostomes and protostomes in figure 1 here: https://onlinelibrary.wiley.com/doi/epdf/10.1111/j.1558-5646.1999.tb03768.x. I recommend that last paper if you want a general sense of how this stuff works.

Because I've now posted several comments, I'll just roll them all up.
For background on the approaches used to identify most recent common ancestors and a high-level look at how animal taxonomy has been inferred, I suggest Lynch 1999.
I think that there are 2 interpretations of this question. If you are interested in just looking up a single MRCA of well-defined clades, such as lobster and human, here are some approaches:
Easy way:

Look at a tree diagram, e.g. this:

Find the tips that correspond to your species of interest (arthropods for lobster, chordata for humans).

Find where they join together in the diagram (the branch labeled "true coelom").

You have your answer, the MRCA is the group of organisms with a true coelom,  coelomates.
A more involved way using a database

Go to this website.
Find the group of species 1 (arthropods, protostomes, etc. for lobster, chordata, deuterostomes etc. for human)
navigate around until you see the group containing the two groups (in this case listed as "bilateria"). In this case you are looking for the bilaterian common ancestor.

another database

Go to this website.
Point and click your way to a view where you see your 2 clades of interest (arthropods, chordates in this case). See figure.

Find where they join (in this case, it is less certain about the existence of a coelomate common ancestor, so it just says "bilaterians").

In R
We might also want a programmatic procedure in case we need to automate it. The following gives a series of commands in R that do this.
Get a newick/nexus tree from somewhere e.g. here of animals (go "control panel --> export, format = "newick tree" --> download).
For this tree, we can approximate an answer as following:
library(ape)
it = read.tree("~/Downloads/zUtVMYn3wJQBoP7E8fygTw_newick.txt")  # autogenerated filename from the website.
getMRCA(it, c("Homo_sapiens", "Drosophila_melanogaster"))  # drosophila is one of only 2 arthropods in this tree, but arthropods all have an arthropod common ancestor
# prints "200", the node index of the MRCA.

You can use that node index to do things like look up a name in the node.label attribute of the tree (it) structure, or to do other operations on the tree.

This is really helpful. I’ll have a good look at this tomorrow and give you a little more detail on what I’m trying to do.

Good answer.
The common ancestor is not known due to the paucity of the fossil record at this point in evolutionary time. A reconstruction can be made using coalescent theory which is what @MaximillanPress describes. The level of uncertainty at that depth of calculation requires significant expertise. I would re-evaluate the comparison with related organisms that are heavily studied.
Lobsters are decapods (Decapoda) and closely related to insects, moreover they are monophyletic with the Order Insecta. A lot is known about a lot of insects, many being model organisms, such as Drosophila melanogaster. Both decapods and insects exist within the arthropods (Arthropoda), along with spiders and the like. Their ancestors were marine arthropods, trilobites are arthropods. The major diversity of the arthropods was generated in the Cambrian Explosion and a lot of work has gone into identifying these fossils.
Traversing between decapods and insects is as far as I'd go, which is likely the late stages of the Cambrian explosion. There are likely to be putative ancestral fossils - not that you can get any molecular sequence from them however. The closest living ancestor on the insect side is Colembola (beach hoppers).
There are serious difficulties in major molecular calculations between decapods and mammals.

I am interested in some of these notes, not being a well-trained systematist myself. Can you say more about the difficulties in molecular clock calculations? Maybe suggest citations for that, or for the beach hopper "living ancestor"?

For the palaeontology of the Cambrian Explosion the author is Simon Conway Morris, e.g. https://www.pnas.org/content/97/9/4426. The arthropod nuclear phylogeny and approximate dates https://www.sciencedirect.com/science/article/pii/S0960982219304865. This includes a guess at the molecular dating (which I'm not convinced by at all). The tree makes the point (which I'd forgotten) that waterfleas (Daphnia) are the more likely a 'living fossils'. Its the lines between the dark blue lines of the Malacostraca (including lobsters) and brown lines of the hexapods (insects).


text:
ERROR: Failed to allocate 18446744073709551596 bytes
I'm using Eagle 2.4.1 for phasing a vcf.gz file. Whenever I give the command it gives me the following error
Phasing samples 1-1
WARNING: Sample 1 (1-indexed) has a het count of 0
ERROR: Failed to allocate 18446744073709551596 bytes

The command to reproduce the error is as follows,
eagle --vcf myprofile.vcf.gz  --geneticMapFile genetic_map_hg38_withX.txt.gz --outPrefix chr2 --chrom 2

Link for the vcf file https://github.com/shafayetrahat/test/blob/main/NA20298_chr2.vcf.gz
Any help will be apreciated.

It seems like you are trying to phase a single sample without a reference panel, that’s also entirely homozygous? You need a reference panel for that and if all the sites are homs, something has gone wrong in the original vcf file.

It looks like a negative number is being stored internally as an unsigned integer, leading to overflow.

@user438383 but the sample is not homozygous. I've checked it. Can the header of the vcf file cause the problem? Also in the eagle documentation, they will check it from server? Anyway can you give link to the ghch38 reference panel server?

@AlexReynolds didn't understand what you mean. Can you speak more elaborately?

@ShafayetRahat post a few lines from the vcf please. I believe Alex means that there may be a bug in the source code somewhere which means Eagle is confused about how much memory is actually being used.

ok I'm adding my vcf in the question

@AlexReynolds added the link for the vcf.gz file

All of your genotypes are 1|1 according to that snippet. That can't be right.

@user438383 apparently, there are big cases for 1/1, or 1|1. But there are many cases of 0/1,1/0 also. I've checked it.


text:
Use of "baseMean" in DeSEQ2
The basemean is described as the "mean of normalized counts of all samples". My question is, how would we interpret differences in basemean between genes, pertaining to reliability and such? For example, if two genes that have similar p-values and l2fchanges, but one has a much higher basemean, would that gene's data be more dependable?
Sorry if this is a trivial question!

I could write a long story about this but actually I will just link this excellent resource from HBCtraining which dissects the individual steps DESeq2 does, starting from normalization over dispersion estimation (this is where the baseMean comes into play) over model fitting and testing. I think it will clarify the role of the baseMean in DESeq2.


text:
Interpreting imputation result from GLIMPSE
I'm following this tutorial of GLIMPSE for learning. I was expecting some extra SNPS coming from the 1000 genome reference in the resulting .vcf file. Though I understand the phasing in the output resulting .vcf, what I don't understand is the imputed result. Apparently the .vcf before and after has the same number of SNPs. I ran the scripts from here. I'm including the vcf files before and after imputation.
Can anyone explain me where the imputed SNPs are in the output v.cf file. If my concept is wrong can anyone tell me what should we expect in result.

Please share the exact command you used.

Just run the scripts from https://github.com/odelaneau/GLIMPSE/tree/master/tutorial

Please paste them into the question to make it easier for myself and others to know exactly what you did. Questions should be as self-contained as possible, a external links break.Thank you.

It seems that you didn't use the following flag:
 --impute-reference-only-variants   Allows imputation at variants only present
                                     in the reference panel. The use of this
                                     option is intended only to allow
                                     imputation at sporadic missing variants.
                                     If the number of missing variants is
                                     non-sporadic, please re-run the genotype
                                     likelihood computation at all reference
                                     variants and avoid using this option,
                                     since data from the reads should be used.
                                     A warning is thrown if reference-only
                                     variants are found.

You were expecting additional SNPs which were present in the reference panel to be imputed into your targets, but you didn't include the flag --impute-reference-only-variants, so only individual missing genotypes were imputed. Individual missing genotypes ./. are always imputed if they are present in both the target and reference, but if the SNP is present in the reference and not the target, it won't be unless you use that flag.

Doesn't answer my question. Can u explain it clearly?

You mean without --impute-reference-only-variants  it will give me only phased data of the sample? Not the imputed snps from reference?

Can you explain your answer more elaborately? I didn't fully understand the answer .

@ShafayetRahat With all due respect, the function of SE isn't to provide a tutorial service for everything you might not totally understand, especially when you say 'I don't understand' rather than being specific about what you don't understand. I feel like my answer provides a clear answer to the question you gave. Having to respond to multiple additional questions in the comments gets a bit tiresome after a while and doesn't;t motivate me to answer other questions you might ask.

sorry about that. Actually, I'm at learning period. So, it may seem that I was asking for a tutorial service. But, in reality, I'm just stuck in the process. I guess the purpose of the whole stack exchange community is for helping with a problem to solve it. It's not like you are bound to answer the question or something like that. You can choose not to answer the question. Anyway, my point is I'm not treating anyone as my personal instructor. I'm very sorry if any of my conversations hurt anyone. And once again sorry for your inconvenience.


text:
install a CRAN package that imports a CRAN package that imports a Bioconductor package
I have a CRAN package say myPack, that imports two CRAN packages, say cranA and cranB. These two packages import Bioconductor packages, say biocA and biocB. When I try installation of my package with install.packages("myPack") this ends up with error. I have to install biocA and biocB that are imported from both cranA and cranB in order to successfully install "myPack".
Does anyone have any idea how to do in order to successfully install myPack without having to install beforehand biocA and biocB from Bioconductor?
P.S. I have seen answers if you add biocViews to the DESCRIPTION file of myPack, but I think this is not sufficient for successful installation with install.packages(). The installation is successful if in addition you set the repositories with setRepositories() and you choose "CRAN" and "BioC software" and then you use install.packages() where the package inside install.packages() contains  biocViews to the DESCRIPTION file. So it seems that biocViews to the DESCRIPTION file does not automatically connect to Bioconductor repository.

I suggest you post this over at support.bioconductor.org.

From my experience there is no good solutions. Many packages I use offer step-to-step guidance about installing CRAN and bioconductor packages separately in a certain order.


text:
How to make MSA for contact prediction?
I'd like to try a contact prediction method (like TripletRes, CCMpred etc) using a custom set of sequences.  There are about 1.5M sequences, all very similar to the canonical sequence; their distribution reflects the occurence of mutations "in the wild" (they're not clustered or anything like that).  How should I make a MSA in this case?
For example in the original TripletRes paper it looks like the MSA comes from DeepMSA - this searches UniClust30 and UniRef90 using hhblits and jackhmmer.  I don't see any way to use DeepMSA to make a MSA from a fully custom set of sequences, and also it's not clear how to combine a custom set with the UniClust/UniRef search results.  I can make an alignment of just my data with jackhmmer, but I don't know if that's what the other tools expect.  For one thing, the alignment is pretty huge compared to the typical inputs (~million vs typical DeepMSA outputs of a few hundred); for another, it is a lot less variable (most sequences differ in just 1-2 aa, vs completely different species).

That's a LOT of sequences! I'm doubtful you'd be able to get TripletRes or CCMPred to run at all with an MSA depth of 1.5 million. I suggest pruning your dataset to remove redundant sequences (sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA). I'm curious if your Neff (normalized number of effective sequences--as defined in the ResTriplet paper you link) is anywhere close to your number of sequences.
Also, I think you've misunderstood the use of DeepMSA. DeepMSA's main use is to find sequences to include in an MSA. I believe it uses much simpler programs to construct the MSA itself. So you'd use DeepMSA if you wanted to do 1 query sequence --> an MSA of several hundred sequences. Since you already have the (very complete) set of sequences, it wouldn't make sense to use DeepMSA. You could instead just use something like Clustal Omega to generate the MSA. There's also a pretty easy-to-use web server but it's limited to 4k sequences.

or have I misunderstood? Are the 1.5M not the same/related protein? If they aren't, I don't see why you'd want to limit your search to the 1.5M only as a larger search would make for a more accurate contact map.

Yes, it's a pretty big dataset - if you're curious, it's all known coronavirus N protein sequences.  I somewhat overestimated initially - it turns out there are only about 77k unique protein sequences in there - taking into account synonymous codons and cases where only a small fragment was sequenced etc.

"I'm curious if your Neff" - You mean Nf, the same thing that is calculated by calcNf in the DeepMSA software package?  I'm still trying to figure that out, the program to calculate it runs out of memory :)

"sequences which don't add much information, i.e. the ones which differ only 1 or 2 AA" - is it definitely true they don't add information? in this case the majority are just 1 or 2 aa changes, but there are a lot of them.  my hope was that this would be more informative than using a much smaller dataset of homologs from different species

You can use sequence weights to investigate how much information you get from a sequence. From the author of DeepMSA, see https://github.com/kad-ecoli/MSAParser/blob/master/calNf_ly.cpp
This should help you downsize your dataset to include the best sequences.

While it's not "definitely true" that a sequence that differs only 1 or 2 AA don't add information, it's likely that the amino acids which change in such a sequence could already change somewhere else in your MSA. Remember that TripletRes and CCMPred look for coevolution in your MSA, so the more variety you have across sequences, the more contacts you can derive from coevolution and the better structure you can model. I strongly recommend using DeepMSA as-is (it should be able to find the sequences in your dataset anyway) to make a good MSA for contact prediction.

Let us continue this discussion in chat.

You might have some luck with Parsnp. From the description:

Parsnp was designed to align the core genome of hundreds to thousands of bacterial genomes within a few minutes to few hours. Input can be both draft assemblies and finished genomes, and output includes variant (SNP) calls, core genome phylogeny and multi-alignments. Parsnp leverages contextual information provided by multi-alignments surrounding SNP sites for filtration/cleaning, in addition to existing tools for recombination detection/filtration and phylogenic reconstruction.

The tldr is that Parsnp is meant for aligning large sets of very closely related microbial sequences. For example, others in my lab have used parsnp to align over a million SARS-CoV-2 genomes. The only potential concern is that the Parsnp workflow also finishes with a non-optional phylogeny estimation step that takes quite long on such large datasets. Someone should really make it so that you can opt out of that step...
Spoiler/disclaimer, as the current maintainer of Parsnp, that someone is probably me :grimace:
The output of Parsnp is in XMFA. You're probably looking for a fasta alignment, so you can use HarvestTools to convert them:
harvesttools -x sars_align/parsnp.xmfa -M sars_align/parsnp.maf

Both parsnp and harvest-tools are available on bioconda, and I'd highly recommend installing them that way.


text:
RNA-Seq Data Heatmap: Is it necessary to do a log2 transformation of RPKM values before doing the Z-score standardisation?
I am making a heatmap using RNA-Seq data in R. The heatmap shows gene expression values (RPKM) in different brain regions. I have the following code:
library(tidyverse)

library(pheatmap)

library(matrixStats)

read_csv("prenatal_heatmap_data.csv") -> all_data

all_data %>%

   column_to_rownames("Brain Region") -> heatmap_data

heatmap_data %>%

   pheatmap()

Which generates the following heatmap:
I want to do a log2 transformation of the RPKM values so that I can do the Z-score standardisation of the data. I have log2-transformed the values using the following code:
heatmap_data %>%

   log2() -> heatmap_data_log2

heatmap_data_log2 %>%

   pheatmap()

However, when I try to create a heatmap using the log2-transformed values I get the following error:

I have looked at the log2-transformed data and I know that the reason I am getting this error is because in the original dataset, some of the RPKM values were 0. When these 0 values get log-transformed they become -Inf:

^ Original Data

^ log2 transformed data
I am not sure how to overcome this issue. I was wondering if it is necessary to do the log2 transformation of the RPKM values before doing the Z-score standardisation? I have seen that it is conventional to do the log2-transformation of the data.
Any advice is appreciated.

Nice questions. Please remember to upvote an accepted answer

Add a pseudocount of one as in the code suggestion I gave you earlier: Best methods for clustering heatmap based on regions where the genes are most highly expressed?
Log of zero is not defined.


text:
Gene names from components in NMF analysis
This question was also asked on Biostars
I am using Rcppml package in R for my NMF analysis, I have a matrix from single-cell analysis. I have cells scores and I wanted to know how to extract genes instead of the components.
i used the following codes.
t1.cpm<- LogNormalize(data = as.matrix(t1@assays$RNA@counts)) # Log converting the data form seurat obj

t1.model <- RcppML::nmf(as.matrix(t1.cpm), 30, verbose = F, seed = 1234)

t1.w <- t1.model$w
t1.d <- t1.model$d
t1.h <- t1.model$h

dim(t1.w)
rownames(t1.w)<- rownames(t1.cpm)
colnames(t1.w)<- paste0("component", 1:30)
rownames(t1.h)<- paste0("component", 1:30)

and my output is as follows
        T_1_TACCTATCAGATGGCA T_1_CCTATTATCAGGCAAG
component1             0.000000000            0.0000000000
component2             0.000000000            0.0023590277
component3             0.001512931            0.0000000000
component4             0.000000000            0.0021574559
component5             0.001175799            0.0020222144
component6             0.000000000            0.0004463346

It would be helpful if you can let me know how to extract Gene names instead of components.
Thanks,
Dave.

Rcppml is a machine learning library.
https://github.com/zdebruine/RcppML
NMF is Non-negative matrix factorization which is used in bioinformatics and AI, which can examine the latent relationships in experimental data sets.

NMF of a gene expression matrix into $k$ components works to decompose the expression matrix into an activity matrix ($H$) and latent factor matrix ($W$).
$W$ describes the $k$ latent factors (often called components) by providing a score for how much each gene contributes to each component. $H$ describes how much each component contributes to the expression of each cell.
To extract the activity and component matrices from your NMF model fit using RcppML, use the following code:
# cells is a Seurat object
model <- RcppML::nmf(cells@assays$RNA@data, k = 5, verbose = F, seed = 1234)
W <- model$w # Components Matrix
H <- model$h # Activity Matrix

# Name genes/components for latent factor matrix
colnames(W) <- paste0("component", 1:k)
rownames(W) <- rownames(cells) # Gene names

# Name components/cells for activity matrix
colnames(H) <- colnames(cells) # Cell IDs
row names(H) <- paste0("component", 1:k)

> head(W)
              component1   component2   component3 component4   component5
MIR1302-2HG 0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
FAM138A     0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
OR4F5       0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00
AL627309.1  2.927886e-07 9.248074e-07 7.961105e-07          0 2.802758e-06
AL627309.3  3.486284e-08 0.000000e+00 1.916799e-07          0 7.938089e-08
AL627309.2  0.000000e+00 0.000000e+00 0.000000e+00          0 0.000000e+00

> head(H)
                 cell_1       cell_2       cell_3       cell_4       cell_5
component1 0.0001370465 0.0002195795 0.0001628898 0.0002239623 0.0001484157
component2 0.0003604048 0.0002239471 0.0005688147 0.0005242678 0.0004821531
component3 0.0003948121 0.0001485325 0.0002885584 0.0002242272 0.0002865587
component4 0.0001391838 0.0002418495 0.0001361317 0.0001382889 0.0001558255
component5 0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000

People will often use the top $N$ genes to describe a component. The following code gets the top $N$ genes for each component and stores them in a list
# This code assumes the above code block has been already run
require(dplyr) #used for %>% piping
N <- 100
top_genes <- list()
for (i in 1:dim(W)[2]) {
  name <- colnames(W)[i]
  top_genes[[name]] <- sort(W[, i], decreasing = T) %>% head(N) %>% names()
}


text:
Steps for trimming aligned peptides using Trimal
I am trimming the aligned peptide sequences and back translating them to codon using Trimal.
Here i am using following cmd's:
trimal -in  OG0008206.fa.aln.out -out OG0008206.fa.aln.out_1st.trim -automated1 -backtrans OG0008206.fa.codon
trimal -out OG0008206.fa.aln.out_2nd.trim -phylip_paml -nogaps

But i am confused about, should i remove the gaps fisrt and then backtransleted the peptides or Backtransleted first and then remove the gaps. Because i am getting different  results (codon) in two attempts. which one i should do fisrt?
Thanks.


text:
'GCLOUD_SDK_PATH' Error
I have run into the following problem
Error: failed to find 'gcloud' binary; set option or environment variable 'GCLOUD_SDK_PATH'

when executing this line of R code:
avfiles_ls( recursive = FALSE,
            namespace = avworkspace_namespace('bdcat-cohort3-stanislawski'),
            name = avworkspace_name('SOLGeneticTopMed')
)

This works when I run it in a notebook in a cloud environment but not when I am running an R script inside a wdl on Terra.
Is there a way to more fully quality where the namespace is?
Thanks.
Liz

Are you running Terra on GCP?

I assume that you are running Terra on Google Cloud Platform (GCP)?
It is weird, because Terra is supposed to run on GCP - that is what its intended purpose - and finding gcloud is automatic, you don't need to configure it.
Three suggestions:

If the cloud is GCP it is appears to could be a bug with the Broad Institutes platform and should be reported directly on their website. They will know immediately, so if they don't respond immediately its a user problem. This is because they will prioritise system bugs, but user errors they may not respond to.
At a guess the error suggests you just need to point to gcloud within .bashrc or .profile, but you don't need to do this on GCP because gcloud is its modus operandi*.
Run WDL directly outside Terra.

The simplest approach is to install miniconda and use it to install gcloud-sdk package and then install Terra from within conda. Assuming Terra is present on "-c bioconda". At this point Terra is forced to see gcloud.
However, it feels like Terra is being installed on a cloud outside GCP and hence the error.

To those who don't know GCP gcloud = google cloud.

Thank you!  Also, the Terra help desk gave me a few lines to add to my docker to do a setup.  Meanwhile, I have found a different way to run my wdl and Rscript so pursuing that for a few days.  Thanks so much for the advice.

Good it helps. Importantly good to know the Broad Institute actively support Terra. Upvotes would be welcome.


text:
What tool RCSB and AlphaFold use to visualize 3D structures?
I would like to implement on a web page a tool to visualize 3D structures like it is for AlphaFold and RCSB.
A couple of examples:
https://www.rcsb.org/3d-view/7KX9
https://alphafold.com/entry/Q5VSL9
As you can see, if you hover over a part of the structure or on a residue on the RNA/amino-acid sequence, that residue is underlined and different details pop up on the screen.
What is the javascript library that they use?
Do you know also where to find the documentation?
Best regards!

Mol*

Different criteria give different rankings.

Mol* (read molstar with trilled rhotic R according to the given IPA) is the newest, is used by the RCSB PDB and can support huge complexes. It is less implemented and has a tricker documentation.
NGL is the former viewer from RCSB PDB. It is good and the switch to Mol∗ was driven by an effort for uniformity with the PDBe. Both are written by Alex Rose. There is more documentation about NGL —in the form of a gallery of codepens and various GitHub repos using it. There is a package for Python notebooks called nglview. It has issues such as no SS predictions (i.e. requires sheets and helix entries to be defined in the file for proper rendering).
3DMol is simpler and is of note as the Python module nglview does not work in colab notebooks, but the Python module py3DMol does.

Great answer! May you apply an edit to your answer and provide links to documentation? I can also do it. In any case, after this addition I will accept the answer.

@Firingam — Totally, made it a community Wiki.

The only caveat about finding the links to the documentation is that the good documentation is generally about the viewers and not the actual JS or TS, for which the GitHub files or issues are the place to look...

There is also the caveat that the use of NGL over Mol* will result in reviewer 2 asking why was the old one chosen over the newer one —been there. :zany_face:


text:
VCF file inspector app with GUI
Does anyone know a good, free-to-use VCF file inspector app (Mac OS, with GUI) where I can browse the content of a VCF file similarly to bamseek*?
I would like to browse SNP annotation split by columns, possibly even filter it.
'* Bamseek does not work for me, fails at open file. It anyways seems to be an abandonware.
Example table view from an expensive commercial software:

I wrote vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html
bcftools view in.bcf |  java -jar dist/vcf2table.jar | more

and http://lindenb.github.io/jvarkit/SwingVcfView.html

java -jar dist/swingvcfview.jar indexed.vcf.gz

OK it's a command line program, not a GUI, but bcftools query and bcftools view can easily be used to to view elements of a vcf/bcf file and are (probably) a lot more flexible and scalable than a GUI. Examples from here:
$ bcftools query -f '%CHROM %POS %REF %ALT\n' file.bcf | head -3
1 13380 C G
1 16071 G A
1 16141 C T

$ bcftools query -f '%CHROM %POS  GTs:[ %GT]\t PLs:[ %PL]\n' file.bcf | head -3
1 10234  GTs: 1/1 1/1    PLs: 28,3,0 29,3,0
1 10291  GTs: ./. 1/1    PLs: 0,0,0 28,3,0
1 14907  GTs: 0/1 0/1    PLs: 8,0,17 26,0,48


text:
Independent Subset of Rectangles Problem for Nonoverlapping local alignments
I'd like to implement a graph-based algorithm for DNA comparison, more especially by solving a problem that can be formulated as a Maximum Weighted Independent Set problem.
I have found the article in reference [1], especially part 2.2.2 "The Independent Subset of Rectangles (IR) Problem", but it is not very explicit on the way to build the graph G given 2 DNA sequence
So I'm looking for any indication, article, or description on the way to construct the conflict graph cited in article reference [1] and [2].
I'm not from the biology world but the computer world.
Thank you very much for your help.
[1] Efficient combinatorial algorithms
for DNA sequence processing: https://www.cs.uic.edu/~dasgupta/resume/publ/papers/book-chapter-dasgupta-kao.pdf
[2] Nonoverlapping local alignments (weighted independent
sets of axis-parallel rectangles) : https://core.ac.uk/download/pdf/82057963.pdf

Consider the below image from a recent review on whole genome alignment:

The lines in the figure are local alignments (sometimes referred to as anchors in literature). These lines can be maintained a number of ways, but often via pairwise genome alignment methods like LastZ, MUMmer, and Last (there are many many more).
Often times, one wants to pull out the "best set" of local alignments. In some cases, algorithms will use the "single copy" heuristic, which constrains that each nucleotide in the query can be aligned to at most one nucleotide from the reference. In this case, you can create a bounding rectangle along each of the lines, and two rectangles are independent if there is no point within either rectangle that is on the same $x$ or $y$ coordinate as any point in the other rectangle. You can use this to set your constraint graph.
However, it is worth noting that two rectangles may overlap just a tiny bit. In this case, it would be ideal to chop the overlapping part off of either rectangle, that way you can select both and not just one. This requires a bit more methodology, but a crude way would be to just look for such cases and for each one, chop the overlap off of the less scoring rectangle.


text:
Is there any function in Biopython to convert a DNA sequence from ambiguous to unambiguous?
I have a project to write in Python that requires to write a function that given  degenerate DNA sequence (for example: KKGTACACCAG) sequence and a molecular weight interval, returns a list of all unambiguous sequences represented by the given sequence (encoding based on IUPAC alphabet: http://www.bioinformatics.org/sms/iupac.html).
Is there any function in Biopython to convert a DNA sequence from degenerate to unambiguous base pairs?
Also, Is it possible to use the molecular weight function from Bio.SeqUtils library for an ambiguous DNA?

from Bio.Seq import Seq
coding_dna = Seq("AAYGANAGYCARAGYAAR")
print(coding_dna.translate())

This gives

NXSQSK

So .translate() will read a degenerate code and provide the unambiguous amino acid, except where it is unable to do so e.g. GAN because this is a split between asp and glu. Thus GAR or GAY is fine, but GAN gives X
You would then use itertools by mapping (map) against the triplet code (import Bio.Data.CodonTable) to obtain all permutations. Doesn't work for X admittedly and stop codons a special case (3 possibilities).
Bio.Data.CodonTable.list_possible_proteins would solve the X problem, but stop codons remain a special case.

I seriously did like the Stackoverflow link in the comments, with the double map funtion ... very slick.

I am missing the SO link the one with :  the double map funtion ... , could you add it to your answer

@pippo1980ONSTRIKE I can't do it now, but I'll look at it early next week. Please check the answer below from acvill

cheers, but R is like ancient greek to me, I am following this question though

M__'s approach is good, but, as explained, some subset of the unambiguous DNA sequences will not be recovered because of the ambiguity in reverse-translation due to the degeneracy of the genetic code. I don't know a python approach, but here's an R function that returns a vector of all the possible unambiguous DNA strings:
make_unambiguous <- function(dna) {
  require(tidyverse)
  require(S4Vectors)
  iupac <- tibble(code = c("A", "C", "G", "T",
                           "R", "Y", "S", "W",
                           "K", "M", "B", "D",
                           "H", "V", "N"),
                  base = c("A", "C", "G", "T",
                           "AG", "CT", "GC", "AT",
                           "GT", "AC", "CGT", "AGT",
                           "ACT", "ACG", "ACGT"))

  tibble(dna) |>
    separate_rows(dna, sep = '(?<=.)(?=.)') |>
    left_join(iupac, by = c("dna" = "code")) |>
    pull(base) |>
    str_split("") |>
    expand.grid(stringsAsFactors = FALSE) |>
    unite(col = sequence, sep = "") |>
    as_vector()
}

make_unambiguous(dna = "AAYGANAGYCARAGYAAR")


text:
Plant signaling pathway interconnections
I'm sorry for having such a juvenile question, but I'm still new in bioinformatics, and I was wondering whether there is any tool able to perform the task I'm dealing with now.
Here is a short description (please let me know if more information should be provided from my side):
I started with 3 genes of Arabidopsis thaliana (let me call them geneA, geneB, geneC). All of them are well annotated and they're known as transcription factors, but each is coming from different signaling pathway.
My aim is to do this:

To find all the signaling pathways related to each of these 3 genes (this is to identify which genes are the most upstream ones in each pathway*).
Afterwards, let's say from all the found signaling pathways of the geneA, I only want to pick those that are sharing genes with the signaling pathways of the geneB and with the signaling pathways of the geneC.

At the end, I would like to get a graphical visualization** of whether and how these signaling pathways are interconnected and what are the particular genes that are shared by those pathways.
*) I cannot use the very popular tool called "String" from string-db.org because it doesn't show the sequence/order of the signaling.
**) I know that KEGG is a database full of graphs of different metabolic patways, but in case of interconnections between signaling pathways in plants, it is either very poor tool, or I'm not that much able to find such a option there to do the job for me (even my colleagues didn't manage finding it).
All of this would be probably more of a literature-mining type of work, but I'm quite sure that, at the end, my task will be to process not only 3 genes, but maybe even more than 10-20 genes. Thus, I would appreciate any more convenient way rather than just digging the info about each gene from the literature and drawing my graphs manually by hand. Automation in this case matters!
Thanks for any help!

Is there any tool that will find signaling pathways that include one or more of my query genes and then graphically visualizes the possible interconnections between those pathways?

i suggest reading more about regulatory network inference tools and related datasets https://www.frontiersin.org/articles/10.3389/fgene.2020.00457/full

@MaximilianPress As far as I understand, such a tool I was looking for does not exist yet. I hope that soon I could find something like a interactive database of signaling pathways interactions. Anyways, thank you for help!

I checked reactome, and that one does not seem to have plants at all. But then google "arabidopsis pathways" returned this https://www.arabidopsis.org/biocyc/? Could it work?


text:
How to split a Seurat cluster in several subclusters?
I've analyzed my scRNA-seq data and have a couple of Seurat clusters that show more than one cell type in each cluster. (for example, cluster 9 shows both NK and CD4 cells)
How can I split a cluster into subclusters, based on the gene expression?

Increase the clustering resolution parameter to generate more (smaller) clusters, see FindClusters in the Seurat docs. Whether or not this will neatly, split your clusters into subclusters depends on your data, but normally one can easily separate CD4 and NK cells from PBMCs. See also the Clustree approach for determining the optimal resolution.
Alternatively/additionally use a cell type annotation approach that does not require clustering, e.g. Seurat V4 permits reference mapping and there are many other methods too

I've done sub-clustering a few times on my Seurat data sets. The approach I take is to subset the clusters that need to be clustered (i.e. using subset), carry out a clustering of only those cells, then transfer the subcluster labels back to the original dataset. Here's some rough code, which will need to be modified for your specific situation and code preferences:
data.sc[["main.cluster"]] <- Idents(data.sc);
clusterToSubset <- "B cells";
data.sub <- subset(data.sc, subset = main.cluster == clusterToSubset);
## ** [Cluster data.sub] **
Idents(data.sc) <- "main.cluster";
## This should work as long as cell order (i.e. column names) isn't changed
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub)) <-
    paste0("Sub_", Idents(data.sub));
data.sc[["merged.cluster"]] <- Idents(data.sc);

This approach will work with multiple sub-clusters as well by changing the label that is pasted onto the sub-cluster idents:
data.sc[["main.cluster"]] <- Idents(data.sc);
clusterToSubset1 <- "B cells";
data.sub.1 <- subset(data.sc, subset = main.cluster == clusterToSubset1);
## ** [Cluster data.sub.1] **
clusterToSubset2 <- "T cells";
data.sub.2 <- subset(data.sc, subset = main.cluster == clusterToSubset2);
## ** [Cluster data.sub.2] **
Idents(data.sc) <- "main.cluster";
## This should work as long as cell order (i.e. column names) isn't changed
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.1)) <-
    paste0("Sub1_", Idents(data.sub.1));
Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.2)) <-
    paste0("Sub2_", Idents(data.sub.2));
data.sc[["merged.cluster"]] <- Idents(data.sc);

Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.1)) gives the same output as Idents(data.sc, cells=colnames(data.sc) %in% colnames(data.sub.2)).
Both are returning the cellIDs of all cells.
The correct code to generate a list of colnames of data.sc is as follows:
Idents(data.sc, cells=colnames(data.sc)[colnames(data.sc) %in% colnames(data.sub.1)]) <- paste0("sub1_", Idents(data.sub.1))

However, you can directly add the new cluster names to a metadata column in the Seurat object using the following code:
## Set up new metadata column to store subcluster IDs later ##
data.sc$subcluster <- Idents(data.sc)

## subclustering steps (from a previous comment)##
data.sub.1 <- subset(data.sc, subset = seurat_clusters == 1)
data.sub.1 <- FindVariableFeatures(data.sub.1)
data.sub.1 <- ScaleData(data.sub.1)
data.sub.1 <- RunPCA(data.sub.1, verbose = FALSE)
data.sub.1 <- RunUMAP(data.sub.1, reduction = "pca", dims = 1:30)
data.sub.1 <- FindNeighbors(data.sub.1, reduction = "pca", dims = 1:30)
data.sub.1 <- FindClusters(data.sub.1, resolution = 0.5)

## merging the subclusster names back to original clusters ##
data.sc$subcluster[colnames(dta.sc)[colnames(data.sc) %in% colnames(data.sub.1)]] <- paste0("sub1_", data.sub.1$seurat_clusters)
```

If those two commands are producing the same output, it means that the subsets (i.e. data.sub.1 and data.sub.2) are identical, i.e. that the subsetting did not work properly.


text:
Evolutionary speaking, what are the constraints on having new mutations?
Evolutionary speaking, some mutations lead to better fitness of an organism and its adaptation to the environment changes, but I was wondering if there are some constraints on having new mutations. For example, if a gene is needed for survival, I expect to see it more resistant to having new mutations. So, my question is what are these constraints, especially for bacteria?

I think the is a very broad question, do you have a more specific problem in your mind? Or example you are puzzled by?

This is an undergraduate/graduate lecture type question and the answer varies: traditional Darwinists believed a reasonable proportion of mutations would result in beneficial adaptation, neutral theory however would state the vast majority of mutations have no effect on phenotype and only an absolute minority result in beneficial adaptation. It depends on which bacteria are the focus of investigation and the selective constraint (I would assume its antibiotics).
Once the bacterial species and selective constraint are described its a straight-forward empirical investigation using molecular mathematics because the mutations driving antibiotic resistance are all known.


text:
How can I get a column with real gene names in my ballgown analysis?
I am doing RNAseq analysis and I am using ballgown procedure for it.
There is an option in R to calculate differentially expressed genes and use FPKM in calculating differential gene expression, which is the following:
results_genes = stattest(bg_filt, feature="gene", covariate="state", getFC=TRUE, meas="FPKM")

When accessing head(results_genes) there is a column named id that contain gene names.
While some genes contain more or less readable names, others look like gene-LOC100037748, gene-LOC100037751, gene-LOC100037760 etc...
Can I create a column with real gene names that correspond to these? Is there a sort of API or something maybe?


text:
An error with annovar
I want to install annovar but in one particular part I always get this error
i1d18@RBGO-Server2:/data2/RNASeq/Angel/annovar$ perl annotate_variation.pl -downdb -webfrom annovar avsnp147 humandb/
NOTICE: The --buildver is set as 'hg18' by default
NOTICE: Web-based checking to see whether ANNOVAR new version is available ... Done
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.gz ... Failed
NOTICE: Downloading annotation database http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.idx.gz ... Failed
WARNING: Some files cannot be downloaded, including http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.gz, http://www.openbioinformatics.org/annovar/download/hg18_avsnp147.txt.idx.gz
fi1d18@RBGO-Server2:/data2/RNASeq/Angel/annovar$

Googling did not return any clue
Any help please?

My first attempt would be to change the version of the genome build to something more recent than hg18.

I don't think hg18 is the issue, according to their documentation hg18 is supported. However, the commands you used are not exactly the same as they showed in the tutorial. Did you try to do exactly what did suggested? If that still does not work, I would suggest to open an issue on their github.

Changing to hg38 solved the problem, thanks


text:
Assigning subcluster idents to original object
I have a scRNA-seq Seurat object I've analyzed, and I noticed that for some of the clusters, there's more than one cell type.
I've created a subset which and run FindClusters again to label the cell more efficiently, and now I want to "paste" the Idents I've assigned in the subcluster to the original object.
How do I do it?

Hi @ShaniS this question does look like it has been answered. Please consider upvoting and/or accepting an answer

Since you haven't shared any data or code so I will show it using the example given below,
First you need to assign a specific cluster to the cells which are having another type like C_new.
cells_df <- data.frame(row.names=paste0("Cell", 1:10, sep = ""),as.vector(paste0("C",1:10,sep="")))

colnames(cells_df) <- "clusters"

cluster_new <- c("Cell6","Cell2")

cells_df$clusters <- as.character(cells_df$clusters)

for(i in 1:length(cluster_new)) {
j <- cluster_new[i]
cells_df[j,] <- "C_new"
}

head(cells_df)
         clusters
 Cell1       C1
 Cell2    C_new
 Cell3       C3
 Cell4       C4
 Cell5       C5
 Cell6    C_new

And if you have already assigned new clusters to the cells then you can skip the above step. Now save the new cluster belongings into the seurat object using the Idents function
Idents(seurat_obj) <- cells_df$clusters

Then you can perform the differential gene expression as follows, for example if you want to find the differentially expressed genes between cluster C7 vs C_new you can do it as follows,
FindMarkers(seurat_obj,ident.1 = "C7", ident.2 = "C_new" )

Assuming parentObj is the original seurat object, and tmp_subset is your subset, you can
new_metadata <- parentObj@meta.data
tmp_subset_annotation <- tmp_subset@meta.data

# new_metadata
           original_annotation
ATTCGGA-1   TCD4
AAATAAA-1   TCD4-naive
...

# tmp_subset_annotation
           new_cellType
ATTCGGA-1   TCD4-naive
AAATAAA-1   TCD4-naive
...

new_metadata$refined_annotation <- sapply(rownames(new_metadata),
     function(ita) ifelse(ita %in% rownames(tmp_subset_annotation),
            tmp_subset_annotation[match(ita,rownames(tmp_subset_annotation)),"new_cellType"],
            new_metadata[match(ita,rownames(new_metadata)),"original_annotation"]))

Basically you iterate each rownames of the original object, and assign the new annotation, if the UMI matches a UMI in the subset, otherwise you keep the original annotation.

Hi, I am getting an error saying
Error in ans[npos] <- rep(no, length.out = len)[npos] :
replacement has length zero
In addition: Warning message:
In rep(no, length.out = len) :

You should add the first few lines of both metadata tables (e.g. head(object@metadata)). Without seeing the data it is difficult to help.

I added a new question with my actual code - https://bioinformatics.stackexchange.com/questions/18365/integrating-subcluster-idents-with-original-object

You should edit your original question rather than create a new thread. Especially because you are asking exactly the same thing.


text:
Mapping NCBI taxID to divergence time?
I want to map NCBI taxID, specifically a pair of Genbank sequences, to the divergence time estimate of their common ancestor. For example, using a coalescence theory approach predominantly used in phylogenetics.
Do we have a database for such correlation?

I want to map NCBI taxID to their evolution time. Do we have a database for such correlation?

It would be best to edit your question to clarify what you mean by "evolution time" since this is not a well-defined term in molecular evolution. Even so, I think we can provide a clear answer to this question.
Estimating branch lengths between divergent taxa is an area of active research—there is no one method that works best in all scenarios. As far as taxonomy is concerned, resources like the NCBI Taxonomy database are under continuous refinement by the community.
Computing estimated divergence times and mapping those times to taxids would be a non-trivial task even for a single clade. Actively maintaining a database of divergence times across the entire taxonomy would be incredibly challenging, and it's doubtful that any such resource exists.

Daniel, thanks for the comment! Do you think it is conceptually not sound if I try to annotate a taxID pair (e.g. taxID A - taxID B) with their divergence time?

@JohnnyTam I do agree with Daniel in that this is an area of active research and divergence times are difficult to calculate even for a single clade.
That said, one source where I've been able to look for publications / estimations of divergence times for genera that I'm interested in has been http://www.timetree.org/. Note that you shouldn't just use the measured statistics like 'Estimated Time' since those are likely to be misleading, but it's a good place to find associated papers.

Frontiers in Ecology and Evolution did a focused and succinct review looking at 'the Tree of Life' in 2017 here. They particularly highlight the importance of integrating a fossil record within molecular trees. The upshot is simply that increasingly tighter integration between palaeontology and molecular dating is the ongoing research frontier.
The best approach in context is to look at relevant literature. The key challenge in my opinion is the underlying population dynamic because this does heavily affect temporal estimates.


text:
How to remedy a DESeq2 collapsing technical replicates error?
Goal: To ensure "the sum of the counts for [my samples] is the same as the counts in the [samples] columns in ddsColl" after collapsing replicates using DESeq2.
Before performing differential expression (DE) analysis on my RNA-seq data stored in a DESeqDataSet object, I need to collapse the read counts observed for each of my technical replicates. This should be simple enough using the DESeq2 manual and vignettes available as reference.
This is how you collapse replicates:
ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
# examine the colData of the collapsed data
colData(ddsColl)
# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

However, when I try this, I get the following error with no apparent resolution reported in the manual or on online forums:
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

How might I circumvent this error and continue on with the analysis?

Input files:
> head(cts)
            KO1_P1_L001.Counts KO1_P1_L002.Counts KO1_P2_L001.Counts KO1_P2_L002.Counts KO1_P3_L001.Counts KO1_P3_L002.Counts KO2_P1_L001.Counts KO2_P1_L002.Counts KO2_P2_L001.Counts KO2_P2_L002.Counts KO2_P3_L001.Counts KO2_P3_L002.Counts WT_P1_L001.Counts WT_P1_L002.Counts WT_P2_L001.Counts WT_P2_L002.Counts WT_P3_L001.Counts WT_P3_L002.Counts
DDX11L1                      1                  0                  1                  1                  2                  1                  0                  1                  2                  1                  3                  1                 7                 2                 1                 1                 0                 0
WASH7P                     124                144                128                151                102                118                 39                 41                103                 87                105                106               125               120                13                26               104               136
MIR6859-1                    8                  6                  4                  5                  4                  6                  2                  2                  7                  2                  7                  6                 6                 6                 2                 2                 2                 9
MIR1302-2HG                  1                  4                  0                  0                  0                  2                  2                  2                  2                  2                  0                  4                 1                 2                 0                 3                 2                 2
MIR1302-2                    0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0
FAM138A                      0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                  0                 0                 0                 0                 0                 0                 0

> print(coldata)
                   condition sample           run
KO1_P1_L001.Counts    longKO     l1  long_P1_L001
KO1_P1_L002.Counts    longKO     l1  long_P1_L002
KO1_P2_L001.Counts    longKO     l2  long_P2_L001
KO1_P2_L002.Counts    longKO     l2  long_P2_L002
KO1_P3_L001.Counts    longKO     l3  long_P3_L001
KO1_P3_L002.Counts    longKO     l3  long_P3_L002
KO2_P1_L001.Counts   shortKO     s1 short_P1_L001
KO2_P1_L002.Counts   shortKO     s1 short_P1_L002
KO2_P2_L001.Counts   shortKO     s2 short_P2_L001
KO2_P2_L002.Counts   shortKO     s2 short_P2_L002
KO2_P3_L001.Counts   shortKO     s3 short_P3_L001
KO2_P3_L002.Counts   shortKO     s3 short_P3_L002
WT_P1_L001.Counts         WT     w1    wt_P1_L001
WT_P1_L002.Counts         WT     w1    wt_P1_L002
WT_P2_L001.Counts         WT     w2    wt_P2_L001
WT_P2_L002.Counts         WT     w2    wt_P2_L002
WT_P3_L001.Counts         WT     w3    wt_P3_L001
WT_P3_L002.Counts         WT     w3    wt_P3_L002

My workflow:
> dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design = ~ condition)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
   condition      sample           run               runsCollapsed
    <factor> <character>   <character>                 <character>
l1    longKO          l1  long_P1_L001   long_P1_L001,long_P1_L002
l2    longKO          l2  long_P2_L001   long_P2_L001,long_P2_L002
l3    longKO          l3  long_P3_L001   long_P3_L001,long_P3_L002
s1   shortKO          s1 short_P1_L001 short_P1_L001,short_P1_L002
s2   shortKO          s2 short_P2_L001 short_P2_L001,short_P2_L002
s3   shortKO          s3 short_P3_L001 short_P3_L001,short_P3_L002
w1        WT          w1    wt_P1_L001       wt_P1_L001,wt_P1_L002
w2        WT          w2    wt_P2_L001       wt_P2_L001,wt_P2_L002
w3        WT          w3    wt_P3_L001       wt_P3_L001,wt_P3_L002

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
logical(0)

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
Error: all(rowSums(counts(dds[, matchFirstLevel])) == counts(ddsColl[,  .... is not TRUE

From the output of dds and ddsColl, it seems the number of columns decreased from 18 to 9, as expected. The only difference I see between my DESeq2 workflow output is the result of matchFirstLevel (see below for expected result). Why might this be for my dataset, but not the example dataset?
Example workflow:
> dds <- makeExampleDESeqDataSet(m=12)
> dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
> dds$run <- paste0("run",1:12)
> ddsColl <- collapseReplicates(dds, dds$sample, dds$run)
> colData(ddsColl)
DataFrame with 9 rows and 4 columns
        condition   sample         run runsCollapsed
         <factor> <factor> <character>   <character>
sample1         B  sample1       run10   run10,run11
sample2         A  sample2        run5          run5
sample3         A  sample3        run6          run6
sample4         A  sample4        run1     run1,run4
sample5         B  sample5        run7          run7
sample6         B  sample6        run9          run9
sample7         A  sample7        run3    run3,run12
sample8         B  sample8        run8          run8
sample9         A  sample9        run2          run2

> matchFirstLevel <- dds$sample == levels(dds$sample)[1]
> matchFirstLevel
 [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE

> stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

Very fair. I wanted to include my input data for reference and reproducibility, but perhaps it makes the content look too daunting/long? I tried to clean it up a little. I think the outputs of the console commands are still necessary because they provide important info for comparison between observed and expected outcomes.

@Gawain It's best to remove lines like colnames(ddsColl) and read.table(...)  since it doesn't affect other people's understanding of the analysis, presumably it's only there for your benefit. Only include lines of code which are exactly necessary for others to understand the analysis you did.

@user438383 I see your reasoning; I have removed those lines from the content :)

sample is a character vector, so levels(dds$sample) is returning null. Presumably you mean to use dds$sample[1] here instead, which will likely solve the problem.

This worked for me! Indeed, in my dataset, sample is a character vector, and in the example dataset, sample is a factor. That would explain why the example code was not working for me. Replacing levels(dds$sample) with dds$sample[1] in the matchFirstLevel variable creation line fixed the problem.


text:
How fo force nextflow to repeat a process until all values in a particular channels are used up BUT a single value from another channel is needed
I am trying to find a hack for getting a process to run until all emissions in a channel are used up. The problem is that the process also takes as input a second channel that only has a single emission. Therefore, the process only runs once, but I need it to repeatedly use the single emission channel for all the emissions in the other channel. I imagine this is not too uncommon a problem, but I have not found a solution.
To be more clear, here is the actual bioinformatics I am trying to do: I have split fastqs into files for multiple samples and emit these files into a channel. I copy this channel because I want to use the files again later. I then build an assembly from the combined files which I index with bwa. I next pass the indexed assembly (index channel) into a new process that also takes as input ALL of the original fastq files (the ones I put into a duplicate channel earlier). The goal is to align each of the fastq files against the genome. However, there are only a few emissions in the index channel so the process doesn't repeat for all the items in my fastq channel.
To aid you all here is a simple example that may help, but take note that my situation doesn't use value channels so is a bit trickier to deal with.
exampleChannel_multiplethings = Channel.from( [1, 3, 5, 7, 9] )
exampleChannel_onething = Channel.from( [1] )

process foo {
  input:
  val x from exampleChannel_multiplethings
  val y from exampleChannel_onething

  output:
  file 'x.txt' into foo_out

  """
  echo $y $x > x.txt
  """
}

This outputs "1 1" into the x.txt file. This is expected behavior, but I need to circumvent it somehow so that the process is repeated for all the items in the exampleChannel_multiplethings channel. Important note, if I was dealing with single values I could probably use one of the combine operators, but in reality I am trying to pass a handful of files for an indexed genome into the process. I haven't found a combine operator that works, but hopefully I have missed something!
Thanks for any assistance
Actual code (overly lengthy to paste here) can be found at: https://github.com/JHarrisonEcoEvo/nextflow_gbs/blob/main/gbs_nextflow

I discovered a solution using a combination of .collect() and each.
If one collects from a path the result is a value channel of paths. Value channels can be used multiple times, whereas queue channels cannot.
ch_collected = Channel.fromPath( "somepath" ).collect()

ch_collected can be used multiple times but if one tries to input it to a process using the usual file..from syntax one will get all the files input to the process at once. Instead, one must use the each option during input. E.g.,
  input:
    each fastq_file from ch_collected

I hope this helps someone.

Thanks for answering your own quesiton :-).


text:
Batch Effects that are Close to Collinear with Main Effect
Two years ago, I started working on a project uses both RNAseq and ATACseq. It's supposed to be a simple Healthy Control (HC) vs disease study. The sample collection was done in 2 phases, and it was clear that there was a batch effect that we could adjust for.
However, recently, we received additional metadata. Plotting a PCA plot showed that there was another, larger batch effect that we didn't account for--location of sample donation. There are 4 different locations with the disease samples being from any of the 4, but ALL HC came from only one of the locations.
I reran the count data through DESeq2 with this design formula: phase + location + disease. It didn't fuss about collinearity like it often likes to do and then pooped out a big list of DEG in the RNAseq.
I could probably run the DEG through GSEA to see if the results match the disease's previously known signatures for biological validity, but what statistical worries should I have about this design matrix? What justification would I need for this statistical asymmetry? Thanks.


text:
imputed dosage values for vcf files
I have one more question related to genotype file. I submitted the job in michigan imputation server.
and my imputed dosage for chr 22 file is:

22 16050435 22:16050435:T:C T C . PASS
AF=0.00098;MAF=0.00098;R2=0.23769;IMPUTED GT:DS:HDS:GP
0|0:0.002:0.001,0.000:0.998,0.002,0.000 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0.001:0.000,0.001:0.999,0.001,0.000
0|0:0:0,0:1,0,0 0|0:0.001:0.001,0.000:0.999,0.001,0.000
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0.001:0.001,0.000:0.999,0.001,0.000
0|0:0.002:0.000,0.002:0.998,0.002,0.000 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0.001:0.000,0.001:0.999,0.001,0.000 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0
0|0:0:0,0:1,0,0 0|0:0.002:0.000,0.002:0.998,0.002,0.000
0|0:0:0,0:1,0,0 0|0:0.001:0.001,0.000:0.999,0.001,0.000
0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0 0|0:0:0,0:1,0,0

Its coming like this.
wherein the test data provided in your predixcan tutorial looks like this:

22_16494187_C_A_b37 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0

Also it states that the imputed dosage should be encoded on a 0-2 scale representing the number imputed or number of the effect alleles the sample posseses. Do anyone know should I convert my vcf files to the desired input genotype file to train the model. I don't understand in my vcf files which column is representing dosage values Thank you.

Does this answer your question? edit dosage file for various columns

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

When i was still working in the field of genotype imputation, a dosage was the expectation (average) of the (posterior) genotype probabilities, coded as GP in the OP's VCF. I once contributed code to bcftools to convert genotype probabilities stored in the GP field to dosages. It looks like that code may now live in a bcftools plugin called "dosage": https://samtools.github.io/bcftools/howtos/plugins.html

Yes, I used bcftools +dosage function to extract the GP field only. thank you.


text:
How to format my file using linux tools?
I have an imputed dosage for chr 22 file. its content is:
#[1]CHROM       [2]POS  [3]REF  [4]ALT  [5]HG00096_HG00096      [6]HG00097_HG00097      [7]HG00099_HG00099      [8]HG00100_HG00100
22      16051249        T       C       0.0     1.0     1.0     0.0    0.0     1.0     0.0     0.0     1.0     0.0

I want to remove [] and remove anything coming after _ from this file as well as instead of separate 22 16051249 T and C i want to write:
22_16051249_T_C_b37

So the final file should look like:
Id      HG00096 HG00097 HG00099 HG00100
22_16051249_T_C_b37     0       1       1       0       1       0       0       1       0

Can anyone please help me out with this.

Is your file tab-separated?

It's better to ask here question regarding specific steps you are struggling with. You seem to want 1, rename header (which frankly I don't know why you would not do in a text editor if you have just one file) 2, merge some columns 3, add there "b37", which I am not even sure where it comes from. I am sure 1 you manage on your own, 2 might be sneakier, however it already is answered many times on stack exchange (if you can't find one it's fine, just post it), 3 you got to explain better. It's not clear at all. Here we like specific questions.

This is what i have tried so far. I tried to remove the brackets using this command and it worked.
sed -e 's/[[^][]*]//g'
And for adding the 22_1605149_T_C I did this:
awk '{ID=$1""$2" "$3""$4" " +b37; print ID}' genotype.chr22_1.txt > genotype.chr22_2.txt
But it just gives me a separate file with these values and not want I want. I tried to join the columns using R but it doesn't generate a separate column. the command used is:
genotype_chr22_1$ID <- paste(genotype_chr22_1$CHROM, genotype_chr22_1$POS,genotype_chr22_1$REF, genotype_chr22_1$ALT, sep = "_")

@RheaBedi I can't really understand the given content of your example. On the first line is the header with 8 columns of which the first 4 will be merged into 1 called Id. On the second line, you have 8 columns, and like the previous, I can understand it will become 22_16051249_T_C_b37. But what about the third line? And how is that in your final file, the first line has 5 columns while the second line has more?

The third line would continue with the values. The first column represent ID and the rest of the columns are sample IDs which represent HG0096. the file would have 2504 sample IDs. and each first row would represent 22_160549_T_C_b37 and their corresponding values 0,1 or 2

Do you mean that the third line is a continuation of the second line? like there's no newline right?

Assuming that the third line is the continuation of the second line. Here's a one (long) liner:
cat <( \
      paste -d '\t' \
         <(echo "Id") \
         <(head -1 test.tsv \
      | cut -f5- \
      | sed -e '$s/\[[[:digit:]]\+\]//g; s/_HG[[:digit:]]\+//g') ) \
    <( \
      paste -d '\t' \
        <(awk 'NR>1 {print $1"_"$2"_"$3"_"$4"_b37"}' test.tsv) \
        <(awk 'NR>1 {print $0}' test.tsv | cut -f5- | sed 's/\.[[:digit:]]\+//g')) \
> output.tsv

Yes, this code worked for my input data.


text:
Mounting a folder from the host to a Docker, while running it as Cromwell runtime
I have a pipeline, which is written in WDL & runs with the help of Cromwell. I use Docker for containerization.
However, my container is dependent on huge (>100GB databases), which cannot be shipped with the image.
Is there an option to mount a directory from the host filesystem to a Docker while running Cromwell workflow? I found this functionality in snakemake, therefore I assume it is common, but I haven't found it in the docs yet.

Not an expert of WDL but looking at Cromwell documentation, there is a block in the configuration of Cromwell called submit-docker:

you'll need to modify the config block
inside backend.providers in your Cromwell configuration. In
particular, this block contains a key called submit-docker

In this example file, on line 420, it is possible to specify the mount option which is -v /path/you/need/to/mount:/path/in/docker/container (usually they are the same) to Docker command line. Setting it like that should work.

Thanks, I modified config and was able to mount directories!