text:
Convert ensemble genes to gene names in my sigle cell signature matrix annotated with seurat in R?
I have single cell data which I want to annotate with the seurat package, so I ran this code in R:
lung.data <- Read10X(data.dir ="C:/Users/Desktop/SINGLE-CELL/DATA",gene.column = 1)

lung <- CreateSeuratObject(counts = lung.data, project = "lung", min.cells = 3, min.features = 200)

And the output of the features of this seurat object is:
    >rownames(lung)
    [1] "ENSG00000000003" "ENSG00000000419"
 [3] "ENSG00000000457" "ENSG00000000460"
 [5] "ENSG00000000938" "ENSG00000000971"
 [7] "ENSG00000001036" "ENSG00000001084"
 [9] "ENSG00000001167" "ENSG00000001460"
[11] "ENSG00000001461" "ENSG00000001497"
[13] "ENSG00000001561" "ENSG00000001617"
[15] "ENSG00000001626" "ENSG00000001629"
[17] "ENSG00000001630" "ENSG00000001631"
[19] "ENSG00000002016" "ENSG00000002079"

How can I change the row names of this Seurat Object from Ensembl to gene name, knowing that not all the ensemble genes have conversion to gene name, in order to obtain the seurat object with the gene names as row names and cell names as column names?

What is a "gene name"? I checked one of your ENSGs and it looks like a human gene, so does that mean you are only working with human sequences? And, if so, do you then mean the HGNC symbol for the gene? Or do you also need to handle other species?

I think BioMart will do the gene ID -> gene name conversion you're after:
https://www.ensembl.org/info/data/biomart/index.html
To use BioMart to convert the gene IDs to gene names, you can use the BioMaRt bioconductor R package. Installation and running instructions can be found here: ensembl.org/info/data/biomart/biomart_r_package.html
For this query, you will need to use the human genes dataset:
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")

Then use your list of gene IDs as the 'Filter', with list of IDs as the values:
filters = 'ensembl_gene_id'

Most human protein-coding genes have an associated HGNC symbol, and these are used in Ensembl. Therefore, you can select the 'gene_name' attribute to retrieve the human gene names, and this will usually correspond to the HGNC gene names, as stored in the Ensembl database.
More information on gene names can be found in the Ensembl documentation: https://ensembl.org/info/genome/genebuild/gene_names.html


text:
PANTHER FUNCTIONAL CLASSIFICATION in R
Does anyone know of a way to programmatically get panther protein class terms and IDs for a gene list, in a way that the classification structure is preserved (PARENT CLASS + CHILD CLASS)?
i.e. For a given gene, I would want to get the list of classes, in away that I can then use the top classification for downstream analyses.
GENEX, PC00244 zinc finger transcription factor > PC00218 DNA-binding transcription factor > PC00264 gene-specific transcriptional regulator

One approach is to get the gene's sequence (DNA/AA) and use the interproscan API  (e.g., https://www.ebi.ac.uk/Tools/common/tools/help/index.html?tool=iprscan5) with the httr2 library in R to POST gene sequence(s) and GET the resulting Panther annotation. Make sure to include Panther in the appl parameter for the post request.
Here's a function that does this for a single sequence.
submit_fasta <- function(
    fasta,
    path_outfolder,
    email,
    stype = "n",
    appl = c("Panther"),
    title = 'stackoverflow'
    ) {
  stopifnot("stype must be 'n' for nucleotide seqs or 'p' for protein seqs" = stype %in% c('n', 'p'))
  url_post <- "https://www.ebi.ac.uk/Tools/services/rest/iprscan5/run"
  chr_seq <- as.character(fasta)
  req <- httr2::request(url_post)
  req <- req |> httr2::req_body_form(
    "email" = email,
    "title" = title,
    "goterms" = "false",
    "pathways" = "false",
    "stype" = stype,
    "appl" = paste0(appl, collapse = ","),
    "sequence" = chr_seq
  )
  req <- req |> httr2::req_headers("Accept" = "text/plain", "Content-type" = "application/x-www-form-urlencoded")
  resp_post <- req |> httr2::req_perform()
  chr_retrieval_code <- resp_post$body |> rawToChar()
  url_get <- "https://www.ebi.ac.uk/Tools/services/rest/iprscan5/result"
  req <- httr2::request(url_get)
  req <- req |> httr2::req_url_path_append(chr_retrieval_code)
  req <- req |> httr2::req_url_path_append("tsv")
  req <- req |> httr2::req_retry(
    max_tries = 15L, backoff = ~ 30L,
    is_transient = function(resp) {httr2::resp_status(resp) == 400L}
  )
  req <- req |> httr2::req_headers(Accept = NULL)
  req |> httr2::req_dry_run(quiet = FALSE)
  resp_get <- req |> httr2::req_perform()
  cat('status code:', resp_get$status_code, '\n')
  result <- resp_get$body |> rawToChar()
  result <- read.table(text = result, sep = '\t')
  return(result)
}

Importantly, the fasta input parameter is expected to be XString object from the Biostrings package. However, the code could easily be adjusted. The sequence eventually needs to be coerced to a string for the POST request.
#fasta <- Biostrings::readDNAStringSet('seq.fa')
#resp <- submit_fasta(fasta = fasta_single, path_outfolder = 'x', email = 'test@gmail.com', title = 'stackoverflow')

edit: set the is_transient param to handle the actual transient error case where the data is not yet available (response code = 400). Also, it's important to note there's a "status" endpoint https://www.ebi.ac.uk/Tools/common/tools/help/index.html?tool=iprscan5#/Status) which is probably much more appropriate than pinging the data endpoint until either success or reaching max tries. However, I hope this answer suffices!

Wow, thank you for your comprehensive answer and welcome to the site.

@M__ Thanks for the welcome reply!

source: Stack Exchange

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

sort: votes


text:
Does AlphaFold actually calculates the fold?
After reading through papers and the code, I'm still not clear if AlphaFold and the likes do actually calculate the 3d structure. My current feeling they are more like information retrieval engines. They find statistical 'similarity' features between the primary structures of the query protein and a database protein (or a family of those) with a known 3d structure, and that known fold is the 'calculated' 3d structure. This would mean that if there is no 'similar' protein with a known 3d fold then the program cannot 'predict' anything.
I realize this could be a stupid layman question but here it goes.

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 you say "calculate", what do you mean? is it "retrieval" when you summarize across a family of similar structures following a sophisticated evolutionary search, comparison, structure-sequence mapping computation, etc.? What's the degree and sophistication of that summary that would mean "calculate" and not "retrieve"?

Your suggested formulation of 'retrieval' is very good. Using molecular dynamics to obtain the 3d structure would be a 'calculation'. My conceptual understanding of AF magic is that it 'calculates' the similarity between the primary structures of the query and a DB target with a known 3d structure ( (this calculation is non-trivial, and seems actually like the crux of the matter). Because similar sequences have similar folds, the 3d structure of the DB target is the answer. It is 'retrieval' because the answer is already known but getting it from the database is hard.

'''the 3d structure of the DB target is the answer''' - this appears to be a fallacy

Alphafold calculates the "folds" throughout a protein, thus its protein structural prediction. It uses convolution-based calculation, which is a modification of co-variance. The reason it doesn't use co-variance directly is complicated, but can result in artefacts.
Thus it looks for "simultaneous" amino acid mutations in protein alignments and believes those represent close spatial proximity - based on previous supervised learning of proteins, comprising an alignment and a known structure.
It has won two CASP competitions (world competition for protein folding prediction). The training is likely deeper both in data sets and ANN/CNN/RNN "layers" (deep learning jargon) than previous protein structural prediction models and hence its success. "Deepmind" the company behind it will be Tensorflow experts, rather than structural biologists.
I don't use it but have seen its structural accuracy first-hand. The accuracy statistic RMSD was impressive.
It is important to note it is built on the back of a long term effort by this community as to calculating protein structural prediction. So it's not "radical technology", they've generally increased the accuracy of established methods. Again they've likely done this with more extensive training.

'''AF calculates the "folds" throughout a protein''' - if this means finding 'piecewise/partial folds' with subsequent 'concatenation' into the final 3d structure, then it would be a 'calculation' indeed.

Thats the gist @biocodz

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


text:
Get a certain gene sequence from bam/vcf and reference
I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I know that vcftools can convert bam to fasta, but I do not see if I can do it only for a certain coordinates (like, for example, chr1:10100500-10200500). I converted the whole bam to fasta and then I tried using bedtools getfasta, but my system ran out of space. Is there a more rational way to do it?

The first approach that comes to mind is to cut the bam file and then convert the already cut bam to fasta:
samtools view -hb chr1:10100500-10200500 > small.bam
samtools fasta small.bam > small.fa

(note that I am not certain about the samtools fasta command, I have never used it and I don't know if it will produce many sequences, one for each haplotype supported by reads; please test this.)

https://www.htslib.org/doc/samtools-fasta.html -- the documentation for samtools fastait should work as @terdon wrote it but just in case you need more information


text:
How to perform liftover from 38 to 37 in R?
I have some gwas summary statistics in GRCh38 that I want to lift to GRCh37. I am trying to liftover in R using this code:
library(tidyverse)
library(magrittr)
library(data.table)
library(rtracklayer)
library(GenomicRanges)

rm(list=ls())

gwas_data <- fread("/gwas_sumstats_allchr.txt")
chain_file <- "/chain_files/hg38ToHg19.over.chain"
chain <- import.chain(chain_file)

# Convert to GRanges object (assuming GENPOS is 1-based)
gwas_ranges <- GRanges(
  seqnames = Rle(gwas_data$CHROM),
  ranges = IRanges(start = gwas_data$GENPOS - 1, end = gwas_data$GENPOS) # Convert to 0-based start
)

gwas_ranges$seqnames <- paste0("chr", gwas_ranges$seqnames)

# Perform the liftover
liftover_output <- liftOver(gwas_ranges, chain)

# Extract successful liftover
liftover_success <- liftover_output$`output`

# Convert the results back to a data frame (convert positions back to 1-based)
result <- data.frame(
  CHROM = as.character(seqnames(liftover_success)),
  GENPOS = start(liftover_success) + 1  # Convert back to 1-based position
)

# Print or save the result
print(head(result))

However this doesn't work and I am getting:
liftover_output <- liftOver(gwas_ranges, chain)
Discarding unchained sequences: 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1, 20, 21, 22, 2, 3, 4, 5, 6, 7, 8, 9

My chain file seems fine on checking it (downloaded from UCSC genome browser), I align my chromosomes to have "chr" like the 38 chain does, and my gwas data is indeed in GRCh38 - what else could be going wrong?
My gwas_ranges looks like:
head(gwas_ranges)
GRanges object with 6 ranges and 1 metadata column:
seqnames      ranges strand |    seqnames
<Rle>   <IRanges>  <Rle> | <character>
[1]       10 10708-10709      * |         chr
[2]       10 10904-10905      * |         chr
[3]       10 10942-10943      * |         chr
[4]       10 12112-12113      * |         chr
[5]       10 13850-13851      * |         chr
[6]       10 14537-14538      * |         chr
-------   seqinfo: 22 sequences from an unspecified genome; no seqlengths

I also tried doing the liftover in python using hail. This ran without error, but the liftover of SNPs on chromosome 10 position 1234 would lift to chromosome 18 12369 . Here's my hail attempt:
import hail as hl
hl.stop()
hl.init()

# Load the reference genomes and add the liftover chain file
rg38 = hl.get_reference("GRCh38")
rg37 = hl.get_reference("GRCh37")
local_chain_file = "/references_grch38_to_grch37.over.chain.gz" #downloaded from Hail's GCP
rg38.add_liftover(local_chain_file, rg37)

# Import GWAS data
gwas = hl.import_table('/gwas_sumstats_allchr.txt',
                       types={'CHROM': 'str', 'GENPOS': 'int'})

# Standardize chromosome names in your GWAS data to match GRCh38
gwas = gwas.annotate(CHROM=hl.if_else(hl.literal('chr').contains(gwas.CHROM), gwas.CHROM, hl.literal('chr') + gwas.CHROM))

# Perform the liftover from GRCh38 to GRCh37
gwas = gwas.annotate(
    locus_grch38=hl.locus(gwas.CHROM, gwas.GENPOS, reference_genome='GRCh38'),
    locus_grch37=hl.liftover(hl.locus(gwas.CHROM, gwas.GENPOS, reference_genome='GRCh38'), 'GRCh37')
)

# Add new columns for CHROM_38 and GENPOS_38
gwas = gwas.annotate(
    CHROM_38=hl.if_else(hl.is_defined(gwas.locus_grch37), gwas.locus_grch37.contig, hl.null('str')),
    GENPOS_38=hl.if_else(hl.is_defined(gwas.locus_grch37), gwas.locus_grch37.position, hl.null('int'))
)

gwas = gwas.drop('locus_grch38', 'locus_grch37')

# Export the results
gwas.export('/gwas_sumstats_allchr_GRCh37.txt')

Here's some of the output I get from the hail run:
Head of gwas_sumstats_GRCh38txt:
CHROM   GENPOS  ID  ALLELE0 ALLELE1 A1FREQ  INFO    N   TEST    BETA    SE  CHISQ   LOG10P  p   MAF
10  10709   10:10709:G:A    G   A   0.0140284   0.420904    56210   ADD 0.143464    0.0384744   13.904  3.71581 0.000192393324862147    0.0140284
10  10905   10:10905:G:A    G   A   0.0352519   0.418554    56210   ADD -0.000990946    0.0246439   0.0016169   0.0141582   0.967925206925404   0.0352519
10  10943   10:10943:G:C    G   C   0.0653494   0.469298    56210   ADD -0.0130032  0.0173743   0.56013 0.342745    0.454208230910946   0.0653494
10  12113   10:12113:AC:A   AC  A   0.0390721   0.431839    56210   ADD 0.0386709   0.0230572   2.81293 1.02915 0.0935082652260488  0.0390721

Head of gwas_sumstats_GRCh37.txt:
CHROM   GENPOS  ID  A1FREQ  INFO    N   TEST    BETA    SE  CHISQ   LOG10P  p   MAF original_locus
18  10905   18:10905:G:A    0.0140284   0.420904    56210   ADD 0.143464    0.0384744   13.904  3.71581 0.000192393324862147    0.0140284   chr10:10709
18  11258   18:11258:G:A    0.0352519   0.418554    56210   ADD -0.000990946    0.0246439   0.0016169   0.0141582   0.967925206925404   0.0352519   chr10:10905
18  11296   18:11296:G:C    0.0653494   0.469298    56210   ADD -0.0130032  0.0173743   0.56013 0.342745    0.454208230910946   0.0653494   chr10:10943
18  12467   18:12467:AC:A   0.0390721   0.431839    56210   ADD 0.0386709   0.0230572   2.81293 1.02915 0.0935082652260488  0.0390721   chr10:12113

The BCFtools/liftover plugin here is exactly aimed at handling GWAS summary statistics and it will do a much better job at dealing with summary statistics for indels


text:
Are there applications of semi-supervised learning to genomics?
Semi-supervised learning is when a portion of the training data set is labeled, and the rest is unlabeled; the unlabeled data is used to improved learning performance / accuracy.
Are there (i.e., have there been) applications of semi-supervised learning to genomics? If so, what are they? If not, is there a reason why not?

Don't mean to be rude, but I literally googled 'semi-supervised learning genomics' and there was loads of relevant hits. Did you try this?

Personal view Semi-supervised learning is basically unsupervised learning with weights, so the "label" or "category" can be quantified rather than approximately clustered as occurs in unsupervised learning.
Unsupervised learning is PCA and tSNE and other multi-variate approaches. The advantages of semi-supervised learning is unlike unsupervised learning the feature in the data that shows greatest proportionate variance can be precisely identified and all features ranked; thus there aren't loads of components where the result could be anywhere.
Sure I've interpreted semi-supervised learning and found it was not useful for any comparative purposes. Basically the same variance Eigen-vector weight observed in one gene in one unique environment, is the strongest Eignen-vector weight for the same gene in a completely different environment. However, these are two independent semi-supervised learning analyses. Thus, it's not measuring comparative difference but relative variance between all genes in that given experiment.
Thus I personally didn't find it useful for my metric, HOWEVER such a result is good for some people where the overall structure of the variance identified by Eigen vectors is very similar  between contrasting environments. Basically it's a comparison to all other genes in the analysis. To get it working you'd need to combine all genes in all environments into a single data set.
That would work, but ensure that the label for each feature (e.g. gene expression) comprised the environment and the gene. What the statistical theoretical validity of such a result would be - I haven't a clue - but it should give some interesting results. To be honest, if I'm not sure about the theoretical validity, I will not be alone.


text:
How long should it take to cluster Uniref?
I have the Uniref90 dataset containing around 76M sequences. I have made modifications to all sequences based on a reduced amino acid alphabet.
I am trying to cluster it again at 90% identity. I am using  mmseq cluster on a 64 core, 500GB RAM machine. It has been more than 12h now.
Starting prefiltering scores calculation (step 1 of 1)
Query db start 1 to 74390789
Target db start 1 to 74390789
[=

How long does it normally take to cluster such a database? I know I can use linclust instead to cluster in linear time, but it might not be what I need in this case.
I am using mmseq cluster as follows:
mmseqs createdb "uniref90-$1.fasta" "uniref90-DB-$1"
mmseqs cluster "uniref90-DB-$1"  "uniref90-DB-$1_clu" tmp/ --cov-mode 2 -c 0.8 --min-seq-id 0.9

Thanks. Hmmm ... when you use htop what do you see when its running?

CPU is used at the maximum and 20% of RAM. Anyway, I have looked into the mmseqs2 in more detail and found that it should take at least a few days.

Thanks and thank you for your question. Thats very clear and answered below.

What does "CPU is used at the maximum" mean? Do you see 100% in the CPU% column? Or do you see 6400%?

6400 %CPU, 12.8 %MEM

When I want to estimate the time taken to do something big, I subset my input data to something that takes a few minutes or less, and calculate forward from there. If there's alignment involved, it's usually a good idea to work it out for multiple dataset sizes, as there could be a linear relationship between input size and time taken, square relationship, or worse.
Have you looked at how long it takes for, say, 1000 sequences, 2000 sequences, and 5000 sequences?

The mmseqs2 code should be specifically compiled and optimised for your system. mmseqs2 optimisation is automated and might be available via Spack* (I'm not sure its actually done for mmseqs2), it is definitely available for EasyBuild. I wouldn't advise attempting a manual optimisation unless you've prior experience.
How much time this will save is more difficult to assess but 50% faster is a conservative assessment. This approach is good practice for a job like this size and its good for the associated report, its standard HPC practice. If you're doing cloud computation (i.e. $$$), HPC optimisation is essential. AWS make serious money from suboptimal HPC practice (in fairness to them they also do a lot of development work to minimise costs).
Finally, if you're an expert in optimisation you can save serious computational time, but this is a major investment of your personal development time and technique. The rationale of how and why is complicated.
*, its very likely available on AWS.

You might simply consider its been running for a few days already so HPC technique isn't needed for this specific job. However, you'll know for the future.

Is the tool multithreaded? The OP mentioned "CPU is used at the maximum" in htop which makes me think that it might be 100% of a single core.

@terdon. Yes it is. Whether it's default - dunno - I was thinking that too. Which is one reason why I requested the htop, but OPs response strongly suggests the parallelisation is default. It's easy to code parallelisation as default (at least in Python), normally its Total CPU - 1. What is a bit unusual I'll admit ... one CPU should be free (standard coding practice) but the OP didn't report that. However, we don't know how the original coder did this. To be honest mmseqs coder is hot, so they'd know the convention. If mmseqs2 is intended for HPC Total CPU - 1 is dropped


text:
Where to access the WGCNA tutorial documents: Horvath lab site down
I am currently using the WGCNA package for some analysis and it seems the Horvath lab site is down. Does anyone know of anywhere else I can access the tutorial documents?

Thank you Peter. Could you also please post the PDF for the FAQs that was originally on the website? They don't seem to be in the dropbox folder. There is a lot of useful tips on there that you and other have referenced to before. Thanks a lot!

Apologies for the web site trouble, we are working with UCLA to restore the web pages. In the meantime, I uploaded the tutorials to Dropbox at https://www.dropbox.com/scl/fo/4vqfiysan6rlurfo2pbnk/h?rlkey=thqg8wlpdn4spu3ihjuc1kmlu&dl=0. I suggest downloading the entire folder, unzipping and opening the file index.html in a browser. You should see a stripped down but usable version of the web page and all tutorial sections should be accessible.

That's great - thanks Peter

Might I suggest using GitHub rather than Dropbox? It's an influential package after all, things like that should be permanently available irrespective of private and academic hosting services.

Yeah, I'm looking into a better and more permanent home for the entire package, GitHub is certainly a candidate...

Would it be possible to get the FAQs section added to the dropbox folder also ?

Any updates @PeterLangfelder?

No, no updates yet :(


text:
(Complete) Stuctural Models of Antibodies: Combining Fc and Fab from pdb to Build Hinge Region? Or what?
I'm a newcomer to the field of antibodies, and I'm seeking guidance on obtaining 3D structures for my molecular simulations. After browsing the PDB databank, I noticed that there are Fab and Fc structures available. My specific interest lies in investigating the structural dynamics in the entire antibody (AB) and on the hinge region.
I have a few questions:

Is it common practice to attempt building the entire antibody structure and hinge regions from PDB files for fragments, or are there recommended alternatives?

Are there comprehensive PDB structures available for different IgG and therapeutic antibodies that include information on the whole structure and hinge regions? Especially for Bevacizumab and Rituximab?

Does it make sense to select single segments, e.g., CH1, VH, VL, for analysis? Would one expect these structures to be reasonably stable?

Additionally, I'd appreciate any insights into common pitfalls and any recommendations you might have for someone new to this field.
Thank you in advance for your time and assistance.
Best regards,
Albert

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

Hi, if you've got multiple questions, it's often better to split them up into separate question posts. However, exploratory posts like this are almost always going to produce better outcomes on discussion websites like reddit

wait for cryomicroscopist to deliver https://www.nature.com/articles/s41467-022-34090-2  Here, we use cryo-EM to image the structure of the human full-length IgM pentamer, r


text:
stacked barplot with error bars
I need to create a barplot with error bars, my data consists of 2 groups called "Genot", with 3 sub-groups "Suberin" for each, associated with a value "Percent".
summary_data <- FY %>%
  group_by(Genot, Suberin) %>%
  summarise(
    Mean_Percent = mean(Percent),
    SD_Percent = sd(Percent))

# plot
ggplot(summary_data, aes(x = Genot, y = Mean_Percent, fill=Suberin)) +
  geom_bar(stat = 'identity', aes(fill = Suberin)) +
  geom_errorbar(data = summary_data,
                aes(x = Genot, ymax = Mean_Percent + SD_Percent, ymin = Mean_Percent - SD_Percent),
                , position = "identity",width = 0.2)

My main problem is that the error bars are placed from 0, but I have to order my data as the following to get the graph in the order I want.
FY$Suberin <- ordered(FY$Suberin,
                    levels = c("Con", "Patchy","Non" ))

I can't figure out a way to describe the positions

I understand why you're trying to do this, but I think that the plot will be kind of hard to look at regardless. I'd suggest making your life simpler and using grouped bars instead.

wait, in your second code block, you are overwriting only that column. You're not actually reordering the rest of the data. Is that intentional?

In my field and for these particular data stacked barplots are what we have to use

In the second code block, I am reordering the full dataset based on the argument "Suberin", this is to obtain my plot with "Non", "Patchy" and "Con" (starting from 0)

(Following from comments)
I am 99% sure that, in the second code block, you are reordering only that one column and the rest of the columns are remaining unordered. That means that that column will be permuted with respect to the rest of the data, which I doubt is what you intend. I'd suggest instead running this, which will reorder the whole dataframe:
# not actually run cause I don't have the data, check it!
FY <- FY[
         order(factor(FY$Suberin,
               levels = c("Con", "Patchy","Non" )
         ),
]

Consider this analogous example:
> test = data.frame(
  "one"=c(1,2,3),
  "two"=c(6,5,4),
  "three"=c("this", "that", "which")
)
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which
> test$three = ordered(test$three, levels=c("that", "this", "which"))
# only column "three" is reordered!!!
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which

# THIS reorders the whole dataframe
> test[order(test$three),]
  one two three
2   2   5  that
1   1   6  this
3   3   4 which

# impose a specific ordering by making it a factor, maybe there's a better way
> test[order(factor(test$three, levels=c("that", "which", "this"))),]
  one two three
2   2   5  that
3   3   4 which
1   1   6  this

I don't know that that impacts your error bar issue. But the underlying data will certainly be affected and possibly made inaccurate by permuting only one column.
Another thing that's a little weird is that you're computing Mean_Percent in two places: first, in your summarise expression, and then again when you actually make the plot with geom_bar(stat="identity"..., which will fill the axes you've set up with newly estimated values based on the relative frequencies of the Suberin column. So you aren't plotting the Mean_Percent values that you are getting from summarise, as far as I can tell.
I'm a little scared of the way that you are computing mean/sd percentages and stacking them, especially not knowing what the data look like. I'd recommend confirming that your Mean_Percent that you compute with summarise adds up to 100%; for the actual bars in the barplot you are letting geom_bar recompute the proportions directly from the data, so I don't think the fact that it apparently adds up to 100% in the plot is trustworthy. To summarize: it'd make sense to confirm that what's in your dataframe matches what's in the plot since they are separately computed.
I don't know exactly that either of these issues is leading to your error bars having the wrong central tendencies, but investigating a couple of these lines may lead you to discover something.
I'm not a ggplot person so I can't diagnose just by looking at your code, and obviously I don't have your data so I can't tinker with it. Someone with more experience there might be a good idea. However, there are some related questions on SO that may be worth reading over.
Note again that none other than Hadley Wickham finds this kind of plot confusing.

Cool answer, @user185980 please consider marking it as accepted. ggplot is a hassle and there's no-one here who'd really do this ... I know it, its just fiddly.

@M__ appreciate the callout, but I'm not sure that my answer really addresses their concern- I didn't solve the problem! I merely have worries about the approach.


text:
Eigen phred score interpretation
I need help. I am not able to find any reference for Eigen_phred_coding prediction tools score from dbNSFP.
I know that it may be scaled in a similar manner to CADD_phred - ranking raw values from highest to smallest; get rank of that value in the queue; phred scale that value with the same formula as sequencing quality score.
But I couldn't find the phred scaling of Eigen values directly, not to mention its suggested cutoff (same as CADD_phred > 15 or different?).
How are the Eigen scores transformed into Eigen_phred?
Is it in same manner as Cadd_phred?
If so, can same cutoff be used in both cases?

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.

@terdon Thank you for that reference. I don't know why, but I always thought that CADD_phred is phred scaled value of ratio rank of value in sorted queue over the total number of used SNV. At least in read me I found out that everything phred is  ratio of rank of the variat over total number of scores in database. Still I am not sure with the cutoff of 15. Anyway I am going to try that.


text:
Generate simulated bam for certain snps
I want to benchmark my DNA sequencing pipeline. In order to do that I want some gold-standard files in every step(.vcf,.bam,.fastq). I want to generate/simulate a bam file of reads for a given set of vcf. How to do that?
The GIAB project has gold standard vcf, bam, and fastq. But they are not perfect for pipeline testing. Their size is quite big. Also, they are real data. My question was how to generate simulated Bam from vcf.

Have you tried making a modified fasta containing the variants and then simulated reads using that. Normally then you'd mix those together with reads simulated against the reference sequence at different ratios.

Actually, I am planning to do that now. The thing is, I couldn't determine which one is the most convenient way. I am planning to do that using Python. Can you write it in the answer and suggest a convenient way to do that?

Hint, you want to make a fasta from a VCF (e.g., https://www.biostars.org/p/360900/) and then use any random read simulator. No need to program anything, that'd be a waste of time.

Sorry, I didn't read your comment earlier properly. Actually, I used many simulators like NEAT, ART, MASON, etc. But, unfortunately, they are too much buggy. If you go to their issues, you will see their problems. That's why I wanted to write my own code in the first place. However, among the simulators, NEAT probably is the least buggy one.


text:
ERROR IN FUNCOTATOR: [SQLITE_IOERR_LOCK] I/O error in the advisory file locking logic (disk I/O error)
I ran the "funcotator" step from the GATK somatic short variant discovery (snvs + indels) pipeline and used the output of filtermutectcalls as input, but I obtain this error:
17:00:18.816 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/gatk/gatk-package-4.3.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:00:19.105 INFO  Funcotator - ------------------------------------------------------------
17:00:19.106 INFO  Funcotator - The Genome Analysis Toolkit (GATK) v4.3.0.0
17:00:19.106 INFO  Funcotator - For support and documentation go to https://software.broadinstitute.org/gatk/
17:00:19.106 INFO  Funcotator - Executing as rita.soares@compute-19 on Linux v5.4.0-148-generic amd64
17:00:19.106 INFO  Funcotator - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_242-8u242-b08-0ubuntu3~18.04-b08
17:00:19.106 INFO  Funcotator - Start Date/Time: November 24, 2023 5:00:18 PM GMT
17:00:19.106 INFO  Funcotator - ------------------------------------------------------------
17:00:19.106 INFO  Funcotator - ------------------------------------------------------------
17:00:19.107 INFO  Funcotator - HTSJDK Version: 3.0.1
17:00:19.107 INFO  Funcotator - Picard Version: 2.27.5
17:00:19.107 INFO  Funcotator - Built for Spark Version: 2.4.5
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:00:19.108 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:00:19.108 INFO  Funcotator - Deflater: IntelDeflater
17:00:19.108 INFO  Funcotator - Inflater: IntelInflater
17:00:19.108 INFO  Funcotator - GCS max retries/reopens: 20
17:00:19.108 INFO  Funcotator - Requester pays: disabled
17:00:19.108 INFO  Funcotator - Initializing engine
17:00:19.577 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FilterMutectCalls/10LC_tumor_S2.filtered.vcf.gz
17:00:19.880 INFO  Funcotator - Done initializing engine
17:00:19.880 INFO  Funcotator - Skipping sequence dictionary validation.
17:00:19.880 INFO  Funcotator - Processing user transcripts/defaults/overrides...
17:00:19.881 INFO  Funcotator - Initializing data sources...
17:00:19.883 INFO  DataSourceUtils - Initializing data sources from directory: /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s
17:00:19.953 INFO  DataSourceUtils - Data sources version: 1.7.2020429s
17:00:19.953 INFO  DataSourceUtils - Data sources source: ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/funcotator_dataSources.v1.7.20200429s.tar.gz
17:00:19.953 INFO  DataSourceUtils - Data sources alternate source: gs://broad-public-datasets/funcotator/funcotator_dataSources.v1.7.20200429s.tar.gz
17:00:20.435 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.annotation.REORDERED.gtf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:20.437 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.pc_transcripts.fa -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.pc_transcripts.fa
17:00:20.513 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:20.584 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Familial_Cancer_Genes.no_dupes.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/familial/hg38/Familial_Cancer_Genes.no_dupes.tsv
17:00:20.711 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/dnaRepairGenes.20180524T145835.csv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dna_repair_genes/hg38/dnaRepairGenes.20180524T145835.csv
17:00:20.779 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_fusion.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_fusion/hg38/cosmic_fusion.tsv
17:00:20.911 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hgnc_download_Nov302017.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/hgnc/hg38/hgnc_download_Nov302017.tsv
17:00:20.967 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/simple_uniprot_Dec012014.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/simple_uniprot/hg38/simple_uniprot_Dec012014.tsv
17:00:21.034 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/CancerGeneCensus_Table_1_full_2012-03-15.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cancer_gene_census/hg38/CancerGeneCensus_Table_1_full_2012-03-15.txt
17:00:21.073 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/achilles_lineage_results.import.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/achilles/hg38/achilles_lineage_results.import.txt
17:00:21.161 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:21.171 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xhgnc_v90_38.hg38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xhgnc/hg38/gencode_xhgnc_v90_38.hg38.tsv
17:00:21.187 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Cosmic.db -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic/hg38/Cosmic.db
17:00:21.341 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:21.357 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xrefseq_v90_38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xrefseq/hg38/gencode_xrefseq_v90_38.tsv
17:00:21.378 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_tissue.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_tissue/hg38/cosmic_tissue.tsv
17:00:21.379 INFO  Funcotator - Finalizing data sources (this step can be long if data sources are cloud-based)...
17:00:21.379 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.annotation.REORDERED.gtf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:21.380 INFO  DataSourceUtils - Setting lookahead cache for data source: Gencode : 100000
17:00:21.553 INFO  FeatureManager - Using codec GencodeGtfCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:21.691 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.pc_transcripts.fa -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.pc_transcripts.fa
17:00:27.190 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.190 INFO  DataSourceUtils - Setting lookahead cache for data source: dbSNP : 100000
17:00:27.412 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.563 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.666 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.814 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Familial_Cancer_Genes.no_dupes.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/familial/hg38/Familial_Cancer_Genes.no_dupes.tsv
17:00:27.837 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/dnaRepairGenes.20180524T145835.csv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dna_repair_genes/hg38/dnaRepairGenes.20180524T145835.csv
17:00:27.853 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_fusion.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_fusion/hg38/cosmic_fusion.tsv
17:00:27.875 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hgnc_download_Nov302017.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/hgnc/hg38/hgnc_download_Nov302017.tsv
17:00:28.142 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/simple_uniprot_Dec012014.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/simple_uniprot/hg38/simple_uniprot_Dec012014.tsv
17:00:28.400 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/CancerGeneCensus_Table_1_full_2012-03-15.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cancer_gene_census/hg38/CancerGeneCensus_Table_1_full_2012-03-15.txt
17:00:28.426 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/achilles_lineage_results.import.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/achilles/hg38/achilles_lineage_results.import.txt
17:00:28.449 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.449 INFO  DataSourceUtils - Setting lookahead cache for data source: ClinVar_VCF : 100000
17:00:28.502 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.702 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.780 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.835 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xhgnc_v90_38.hg38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xhgnc/hg38/gencode_xhgnc_v90_38.hg38.tsv
17:00:29.651 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Cosmic.db -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic/hg38/Cosmic.db
17:00:29.807 INFO  DataSourceUtils - Setting lookahead cache for data source: Oreganno : 100000
17:00:29.828 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:29.871 INFO  FeatureManager - Using codec XsvLocatableTableCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.config
17:00:30.155 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:30.158 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
WARNING 2023-11-24 17:00:30     AsciiLineReader Creating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream
17:00:30.168 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xrefseq_v90_38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xrefseq/hg38/gencode_xrefseq_v90_38.tsv
17:00:30.284 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_tissue.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_tissue/hg38/cosmic_tissue.tsv
17:00:30.398 INFO  Funcotator - Initializing Funcotator Engine...
17:00:30.400 INFO  FuncotatorEngine - Using given VCF and Reference.  No conversion required.
17:00:30.402 INFO  Funcotator - Creating a MAF file for output: file:/mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FUNCOTATOR/10LC_tumor_S2.FiltPASS.maf
17:00:30.420 INFO  ProgressMeter - Starting traversal
17:00:30.420 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
17:00:42.145 INFO  ProgressMeter -         chr1:1129652              0.2                  1000           5117.3
17:00:52.621 INFO  ProgressMeter -         chr1:1544800              0.4                  2000           5405.4
17:01:03.582 INFO  ProgressMeter -         chr1:2088265              0.6                  3000           5427.9
17:01:23.573 INFO  ProgressMeter -         chr1:2868977              0.9                  5000           5644.2
17:01:45.697 INFO  ProgressMeter -         chr1:3708722              1.3                  7000           5579.4
17:01:47.516 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000581512.2 for variant: chr1:3783110-3783112(TAA* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000581512.2  Variant: [VC Unknown @ chr1:3783110-3783112 Q. of type=INDEL alleles=[TAA*, T] attr={AS_FilterStatus=[weak_evidence, contamination], AS_SB_TABLE=[0, 0|0, 0], DP=1, ECNT=2, GERMQ=1, MBQ=[0, 31], MFRL=[0, 259], MMQ=[60, 60], MPOS=66, POPAF=2.21, ROQ=93, RPA=[15, 13], RU=A, STR=true, STRQ=93, TLOD=4.20} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:PGT:PID:PS:SB     0|1:0,1:0.667:1:0,0:0,1:0,1:0|1:3783099_C_T:3783099:0,0,1,0 filters=haplotype,germline,contamination,weak_evidence

... (CODE OMITTED BECAUSE OF CHARACTERS LIMINATION OF THE BODY OF THE QUESTION)

17:06:03.884 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000375071.4 for variant: chr1:20589096-20589097(TC* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000375071.4  Variant: [VC Unknown @ chr1:20589096-20589097 Q. of type=INDEL alleles=[TC*, T] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[108, 95|79, 72], DP=369, ECNT=2, GERMQ=1, MBQ=[20, 20], MFRL=[119, 123], MMQ=[60, 60], MPOS=26, PON=true, POPAF=0.256, ROQ=93, RPA=[3, 2], RU=C, STR=true, STRQ=93, TLOD=206.22} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB   0/1:203,151:0.422:354:57,38:55,42:114,83:108,95,79,72 filters=panel_of_normals,germline
17:06:03.885 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000375071.4 for problem variant: chr1:20589096-20589097(TC* -> T)
17:06:20.587 INFO  ProgressMeter -        chr1:21789645              5.8                 38000           6511.2
17:06:38.318 INFO  ProgressMeter -        chr1:22910935              6.1                 40000           6523.5
17:06:55.667 INFO  ProgressMeter -        chr1:24196391              6.4                 42000           6541.3
17:07:12.288 INFO  ProgressMeter -        chr1:25333744              6.7                 44000           6569.3
17:07:29.628 INFO  ProgressMeter -        chr1:26561968              7.0                 46000           6583.9
17:08:54.624 INFO  ProgressMeter -        chr1:27258140              8.4                 47000           5593.0
17:09:11.523 INFO  ProgressMeter -        chr1:28509457              8.7                 49000           5641.9
17:09:14.750 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000294409.2 for variant: chr1:28714873-28714874(TA* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000294409.2  Variant: [VC Unknown @ chr1:28714873-28714874 Q. of type=INDEL alleles=[TA*, T] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[1, 8|0, 5], DP=14, ECNT=1, GERMQ=1, MBQ=[29, 29], MFRL=[277, 190], MMQ=[60, 60], MPOS=19, PON=true, POPAF=1.96, ROQ=93, RPA=[13, 12], RU=A, STR=true, STRQ=1, TLOD=6.90} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB   0/1:9,5:0.377:14:1,1:4,1:8,5:1,8,0,5 filters=panel_of_normals,germline,slippage
17:09:14.750 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000294409.2 for problem variant: chr1:28714873-28714874(TA* -> T)
17:09:28.156 INFO  ProgressMeter -        chr1:29742106              9.0                 51000           5690.5
17:09:44.363 INFO  ProgressMeter -        chr1:31000060              9.2                 53000           5740.7
17:10:02.660 INFO  ProgressMeter -        chr1:32309371              9.5                 55000           5766.8
17:10:21.404 INFO  ProgressMeter -        chr1:33875394              9.8                 57000           5787.0
17:10:39.604 INFO  ProgressMeter -        chr1:35564548             10.2                 59000           5811.1
17:11:13.558 INFO  ProgressMeter -        chr1:36959618             10.7                 61000           5690.8
17:11:30.930 INFO  ProgressMeter -        chr1:38301837             11.0                 63000           5722.9
17:11:48.409 INFO  ProgressMeter -        chr1:39672835             11.3                 65000           5752.3
17:12:04.997 INFO  ProgressMeter -        chr1:40888483             11.6                 67000           5787.7
17:12:22.171 INFO  ProgressMeter -        chr1:42389559             11.9                 69000           5816.6
17:13:07.894 INFO  ProgressMeter -        chr1:43856329             12.6                 71000           5624.0
17:13:25.958 INFO  ProgressMeter -        chr1:45154774             12.9                 73000           5647.7
17:13:43.676 INFO  ProgressMeter -        chr1:46612646             13.2                 75000           5672.8
17:14:00.829 INFO  ProgressMeter -        chr1:48197446             13.5                 77000           5700.8
17:14:19.050 INFO  ProgressMeter -        chr1:50725644             13.8                 79000           5720.3
17:14:29.607 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000361556.9 for variant: chr1:51729808-51729817(GGGGGGTGAA* -> G): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000361556.9  Variant: [VC Unknown @ chr1:51729808-51729817 Q. of type=INDEL alleles=[GGGGGGTGAA*, G] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[44, 1|46, 0], DP=110, ECNT=1, GERMQ=1, MBQ=[37, 37], MFRL=[163, 167], MMQ=[60, 60], MPOS=35, POPAF=1.89, ROQ=93, RPA=[2, 1], RU=GGGGGTGAA, STR=true, STRQ=93, TLOD=175.52} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:45,46:0.514:91:23,24:12,22:44,46:44,1,46,0 filters=germline
17:14:29.607 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000361556.9 for problem variant: chr1:51729808-51729817(GGGGGGTGAA* -> G)
17:14:29.715 INFO  ProgressMeter -        chr1:51744810             14.0                 80000           5719.1
17:14:47.317 INFO  ProgressMeter -        chr1:53112621             14.3                 82000           5741.6
17:15:03.760 INFO  ProgressMeter -        chr1:54313159             14.6                 84000           5770.9
17:15:20.710 INFO  ProgressMeter -        chr1:55702719             14.8                 86000           5795.9
17:15:38.363 INFO  ProgressMeter -        chr1:57618601             15.1                 88000           5815.3
17:15:55.381 INFO  ProgressMeter -        chr1:59379903             15.4                 90000           5838.1
17:16:13.770 INFO  ProgressMeter -        chr1:61650616             15.7                 92000           5851.5
17:16:29.978 INFO  ProgressMeter -        chr1:63193764             16.0                 94000           5877.7
17:17:46.249 INFO  ProgressMeter -        chr1:64973681             17.3                 96000           5560.8
17:18:29.297 INFO  ProgressMeter -        chr1:65743799             18.0                 97000           5394.5
17:18:44.743 INFO  VcfFuncotationFactory - ClinVar_VCF 20180429_hg38 cache hits/total: 0/138
17:18:44.746 INFO  VcfFuncotationFactory - dbSNP 9606_b151 cache hits/total: 0/63294
17:18:44.748 INFO  Funcotator - Shutting down engine
[November 24, 2023 5:18:44 PM GMT] org.broadinstitute.hellbender.tools.funcotator.Funcotator done. Elapsed time: 18.43 minutes.
Runtime.totalMemory()=2054684672

org.broadinstitute.hellbender.exceptions.GATKException: Unable to query the database for geneName: LINC01702
        at org.broadinstitute.hellbender.tools.funcotator.dataSources.cosmic.CosmicFuncotationFactory.createFuncotationsOnVariant(CosmicFuncotationFactory.java:320)
        at org.broadinstitute.hellbender.tools.funcotator.DataSourceFuncotationFactory.determineFuncotations(DataSourceFuncotationFactory.java:245)
        at org.broadinstitute.hellbender.tools.funcotator.DataSourceFuncotationFactory.createFuncotations(DataSourceFuncotationFactory.java:211)
        at org.broadinstitute.hellbender.tools.funcotator.FuncotatorEngine.createFuncotationMapForVariant(FuncotatorEngine.java:178)
        at org.broadinstitute.hellbender.tools.funcotator.Funcotator.enqueueAndHandleVariant(Funcotator.java:924)
        at org.broadinstitute.hellbender.tools.funcotator.Funcotator.apply(Funcotator.java:878)
        at org.broadinstitute.hellbender.engine.VariantWalker.lambda$traverse$0(VariantWalker.java:104)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.accept(ForEachOps.java:183)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.stream.ReferencePipeline$2$1.accept(ReferencePipeline.java:175)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.Iterator.forEachRemaining(Iterator.java:116)
        at java.util.Spliterators$IteratorSpliterator.forEachRemaining(Spliterators.java:1801)
        at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:482)
        at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:472)
        at java.util.stream.ForEachOps$ForEachOp.evaluateSequential(ForEachOps.java:150)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.evaluateSequential(ForEachOps.java:173)
        at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:234)
        at java.util.stream.ReferencePipeline.forEach(ReferencePipeline.java:485)
        at org.broadinstitute.hellbender.engine.VariantWalker.traverse(VariantWalker.java:102)
        at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1095)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)
Caused by: org.sqlite.SQLiteException: [SQLITE_IOERR_LOCK]  I/O error in the advisory file locking logic (disk I/O error)
        at org.sqlite.core.DB.newSQLException(DB.java:1030)
        at org.sqlite.core.DB.newSQLException(DB.java:1042)
        at org.sqlite.core.DB.execute(DB.java:884)
        at org.sqlite.core.CoreStatement.exec(CoreStatement.java:80)
        at org.sqlite.jdbc3.JDBC3Statement.executeQuery(JDBC3Statement.java:68)
        at org.broadinstitute.hellbender.tools.funcotator.dataSources.cosmic.CosmicFuncotationFactory.createFuncotationsOnVariant(CosmicFuncotationFactory.java:288)
        ... 26 more
Using GATK jar /gatk/gatk-package-4.3.0.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /gatk/gatk-package-4.3.0.0-local.jar Funcotator -R /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/bwa_prefiles/GRCh38.primary_assembly.genome.fa -V /mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FilterMutectCalls/10LC_tumor_S2.filtered.vcf.gz -O /mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FUNCOTATOR/10LC_tumor_S2.FiltPASS.maf --output-file-format MAF --ref-version hg38 --data-sources-path /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s --disable-sequence-dictionary-validation
srun: error: compute-19: task 0: Exited with exit code 3
srun: launch/slurm: _step_signal: Terminating StepId=2949215.0

Are you sure there's a problem? The error appears after the Funcotator done. Elapsed time: 18.43 minutes. message, which suggests it had actually finished. Then, the error itself (Unable to query the database for gene... and Caused by: [...]  I/O error in the advisory file locking logic (disk I/O error)) suggests a problem with reading some DB that may well be transient. Did you try running it again?

Yes i did runned again but i got the same error and i only obtain a vcf file and not a maf file as output of Funcotator


text:
Energy error margin for protein structures
I would like to know the acceptable error margin for the lowest free energy of a protein. In other words, say, both Group A and B have determined the lowest free energy of the same protein in their respective labs. What is an appropriate error margin so that I can compare these two results?
I just came across this article.
Now let me quote the following from the article.

Note that we generally require predictions to be within 1 kcal/mol to
be useful, since 1.4 kcal/mol equates to a 10-fold difference in
binding affinity.

Can we use this 1 kcal/mol as the acceptable error margin?
Edited:
from *(lehninger principles of biochemistry) and Is Kd (the dissociation constant) for a given protein-ligand pair directly correlated to binding affinity (kcal/mol, Rosetta Energy Units, etc.)?

"... both Group A and B have determined the lowest free energy of the same protein in their respective labs .... " do you mean theoretically or empirically ... in both cases I would like to know more about standard ways to achieve that nowadays

talking about "The maximal and current accuracy of rigorous protein-ligand binding free energy calculations" https://www.nature.com/articles/s42004-023-01019-9 ---> ...As our experimental reproducibility survey indicated that experimental RMSE is on average ~1 kcal mol−1, it would be an extraordinary challenge to ever have an FEP method that achieves an error truly statistically indistinguishable from the experimental error o.... So again is the error margin sought a theoretical or empirical one ??

what kind of chemical interaction gives rise to a change of 1.4 kcal/mol ??

@pippo1980, I had empirical determination in mind.

@pippo1980, I do not know what kind of chemical interaction gives such rise. Sorry.

dont worry , just curious about your question , but not reaaly knowledgable about it , you could try https://mattermodeling.stackexchange.com/ too : Matter Modeling Stack Exchange is a question and answer site for Matter Modelers: computational chemists, material scientists, particle physicists, data scientists, and anyone else who uses computational methods to study molecules and materials


text:
issue with rna seq analysis
I am working on RNA seq analysis and I would like to know the following things:
Current Methods:

downloaded genome fasta file for non-coding rna from ensembl and got the gtf file for hg38
performed hist2 and got 17% alignment for my sample against ncrna genome
sorted the file and then ran feature counts with the hg38 gtf file but got zero counts

Question:

what could be the issue that is causing zero counts?

My bam file after sorting has the SRR ID at first whereas the GTF file has chromosome no: . Could this be the issue?

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.

Did you spot-check where the alignments are?

nope. i guess there was an issues with the refernce genomecompatibility.

Several issues:

There is no non-coding RNA genome. There are non-coding genes but the genome fasta files from Ensembl are always the entire genome and this is what you should align against. Probably you downloaded a transcriptome fata file.

If you aligned against a transcriptome fasta file (like ncRNAs) then it makes no sense to use a GTF file since GTF files have genome coordinates but you mapped against a selection of a transcriptome fasta, hence you have (specific) transcriptome coordinates. GTF and these coordinates are not compatible.

Solution: Align your data against the entire genome (hg38), then use the GTF to get counts for all annotated genes, then do subsetting to specific genes if that is what you need.

Thank you so much. I guess i have been using transcriptome fasta file. will change that .


text:
Generating force-fields for ANY residue
Is this possible to generate force fields for any residue in my .pdb file? I am trying to generate conformations for the proteins with arbitrary residues binded to it from .pdb files using OpenMM, but there are always errors like 'missing residue VAL, missing residue HEM, etc'. I tried to use pdbfixer from openmm library, skip\replace missing\nonstandard residues using GROMACS and generate .top .gro, but even after that it shows me errors or generates nonsence if I replace residues. I read some articles and looks like generating force fields is a complicated task. But still wondering is there any open-source framework\library\neural network that can generate such force fields?
For example I follow this guide
http://docs.openmm.org/latest/userguide/application/02_running_sims.html#a-first-example
Whenever I try to run the following code
from openmm.app import *
from openmm import *
from openmm.unit import *
from sys import stdout

pdb = PDBFile('input.pdb')
forcefield = ForceField('amber14-all.xml', 'amber14/tip3pfb.xml')
system = forcefield.createSystem(pdb.topology, nonbondedMethod=PME,
        nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinMiddleIntegrator(300*kelvin, 1/picosecond, 0.004*picoseconds)
simulation = Simulation(pdb.topology, system, integrator)
simulation.context.setPositions(pdb.positions)
simulation.minimizeEnergy()
simulation.reporters.append(PDBReporter('output.pdb', 1000))
simulation.reporters.append(StateDataReporter(stdout, 1000, step=True,
        potentialEnergy=True, temperature=True))
simulation.step(10000)

It shows
No template found for residue (number) (residue_name)

for every .pdb file that I am trying to run and every ForceField parameter permutation.
As an example I tried to run simulations on this protein
https://www.rcsb.org/structure/1a3n

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 you say "force field", can you be more specific? it seems like a force field is a step in your workflow rather than your goal. "generate conformations for the proteins with arbitrary residues binded to it": what does this mean? What exactly are you modeling?  what exactly are you doing to your pdb file? can you give examples?

Do you mean "electron density clouds"?


text:
What are the applications of Tries(data structure) of an ordered sequence of strings in bioinformatics?
This question was also asked on reddit
Tries are a data structure that can be used to efficiently store and search for strings.
Tries created from an ordered sequence of strings  differ from the regular Tries in the following way:
If there is a string of a certain character $w$ at a character position $x$ and is between two other strings that have both the same character $p$ where $p \neq
 w$ at position $x$(and all their characters before position $x$ are the same in all three strings) then the Trie must have an edge for the first occurrence of $p$, an edge for $w$ and an edge for the second occurrence of $p$.
I have recently written a paper regarding a novel data structure that is related to tries and I wish to connect it to some applications in bioinformatics since the paper itself presents only one application related to string algorithms.
Since bioinformatics makes some use of string algorithms and since the structure is related to Tries I wish to ask you about the possible applications of Tries so I may find some application to my related data structure.
What are the applications of Tries(data structure) of an ordered
sequence of strings in bioinformatics?

This seems a bit like an answer searching for a question. Discussion-based questions don't work so well here, and are better asked on discussion sites like reddit. Are you able to clarify your particular problem and relate it to a specific bioinformatics problem?

(Should be a comment not an answer, but I weary of SE comment formatting- feel free to express disapproval in whichever way)
When in doubt, google it and look at wikipedia:

Tries are used in Bioinformatics, notably in sequence alignment software applications such as BLAST, which indexes all the different substring of length k (called k-mers) of a text by storing the positions of their occurrences in a compressed trie sequence databases.[25]


text:
bwa mem hangs after a few thousand reads
I am trying to align a bunch of paired sample fastq files using bwa mem.
My original command was:
bwa mem -t 8 hg38.fa sample_read1.fq.gz sample_read2.fq.gz > sample_paired.sam

I am running this on a HPC cluster.
These files have approx. 25 million reads, so I initially anticipated that they might take a little bit of time; however I noticed that the programs seemed to hang - the sam file that I am redirecting the output to never became populated with any entries other than the header comments, even after hours of waiting.
Noticing this, I decided to first test whether the most bare bones bwa mem functionality was working, and so I took just the first 1000 lines off one of the fastq files and tested whether it was working with:
bwa mem hg38.fa test_1000lines.fq.gz > test_output.sam

So this works. This creates a proper sam file almost instantaneously. However, I kept testing this behavior by increasing the number of lines I took off the fastq file from 1000 to 2000,3000,etc., and eventually realized that any fastq files I input with more than approximately 8400-8500 lines will make the program hang.
If I run the following with a fastq file containing 8400 lines, the program will run and the sam file will be created in a few seconds:
bwa mem hg38.fa test_8400lines.fq.gz > test_output8400.sam

and this is the stdout output:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 2100 sequences (313733 bp)...
[M::mem_process_seqs] Processed 2100 reads in 1.051 CPU sec, 1.054 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem supporting_files/hg38.fa test_8-4k.fq.gz
[main] Real time: 4.231 sec; CPU: 4.178 sec

However, if I run the following with just 100 more lines in the fastq file:
bwa mem hg38.fa test_8500lines.fq.gz > test_output8500.sam

The command will hang (at least, that is what it appears like it's doing) and all I will see on the terminal is:
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 2125 sequences (317468 bp)...

followed by a blinking cursor.
I've monitored the CPU and MEM usage with the "top" command while I was conducting these tests, and while the memory never goes above something like 15%, the CPU usage always goes to 100% when this hanging behavior happens. When the .sam files were created properly and the program didn't hang, the CPU usage was (for the very short amount of time it took for the program to run) like 60-70%.
Whether the CPU usage going to 100% is causing the hanging, or the hanging it causing the CPU usage to 100%, I'm not sure.
I'm completely out of ideas what the issue is here. Any ideas?
EDIT: Although I originally suspected this to be an issue with my fastq file, I observed the stalling behavior with a different fastq file from a completely different source and ruled that out. It turns out the problem is something related to the reference index, since when I use a different set of reference indices (the hs38DH set) instead of the hg38 set that I had been using, the problem disappears.

The problem with this question is there's no diagnosis. top doesn't diagnose the cluster - it diagnoses the head node - where the job usually shouldn't be running by convention (cluster depending). RAM allocation between the head node and the worker nodes isn't always the same thing.  The correct command is htop. Whats the output of that? -t is multi-threading. htop will describe the CPU and RAM allocation across every CPU in the cluster - not just one node. Are you are running a queue system and a inhouse cluster (via jobs) or a cloud and partitions with/without job submission?

Does this always happen or is it specific to these particular reads only?

@terdon So what I've been able to confirm since I posted this question, is that there are multiple locations in this fastq file where this hanging behavior occurs (for example, if I skip over the first 50k lines and take out a chunk of 10k lines and run bwa mem, then I will once again observe stalling behavior from those 10k lines). It looks like this problem between line 8400 and 8500 just happens to be the first occurrence of the problem. I can also confirm that my other sample fastq files also exhibit this behavior, although at different locations in the file.

Please edit that into your question. Also tell us what species you are sequencing, how much of it, if your target includes many repetitive or ow complexity regions etc. It really sounds like you're hitting some specific issue but I cannot be sure and maybe if you add more detail someone more knowledgeable than I can figure it out. Even better, try and get the smallest file that reproduces the issue and upload it somewhere we can see it (assuming you are allowed to).

Cross-posted here: bwa mem hangs after a few thousand reads

@terdon So it turns out that if I use a different set of reference indices (for ex. the hs38DH set) as opposed to the hg38 set which I built using "bwa index", then the problem disappears.
This solves my immediate problem, although I'm curious if you have any insights into why the original hg38 build indices might have caused bwa mem to stall? I have the log from the index build and it does not have any complaints or errors. As far as I can tell, it worked fine.

I can't know without knowing the exact commands you ran, what kind of preprocessing you did to the genome, which minor version of hg38 etc etc.You could try asking a new question with these details but it will be hard for us to debug without access to the actual data.

I mean it was literally just these 3 lines:

wget -P reference/ https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

gunzip reference/hg38.fa.gz

bwa index reference/hg38.fa

did you md5sum the downloaded references? are you sure the download was complete and not corrupt?


text:
Kallisto RNA-seq errors because of incompatible indices when running <kallisto quant> using downloaded mouse index
I'm trying to run RNA-seq on single-end mouse RNA data (50bp reads) but have very limited RAM (<4 gb) - I don't think I really have a way to work around this memory restriction, and wanted to try Kallisto's pseudoalignment (I conda installed it) as indexing and/or aligning with other tools like STAR and HISAT2 have resulted in the processes being killed. Building the index with kallisto index still results in the process being killed, so I downloaded the index files from https://github.com/pachterlab/kallisto-transcriptome-indices/releases/tag/v1 instead. I ran kallisto quant -i index.idx -o output --single -l 180 -s 20 sample1.fastq (I don't entirely know what the mean and sd arguments should be for my data), but this resulted in the following:
[quant] fragment length distribution is truncated gaussian with mean = 180, sd = 20
Error: incompatible indices. Found version 13, expected version 10
Rerun with index to regenerate

I've seen that other people have come across this error before when quantifying with kallisto, but it seems that they are building the index themselves while I was trying to use the pre-generated index from the unpacked mouse_index_standard.tar.xz file downloaded from the link above.
How can I get around this error? And in general, are there any suggestions on how to get raw/normalized count data from RNA-seq data with such memory limitations? I'm rather new to this.

If your finances at all allow, get more RAM. You're not soing much with such limited resources in the long run.

It means that the version of the tool used for indexing and your current version are different and not compatible. Either downgrade your version to the version used for indexing, or build a new index with your current version. As by the GitHub repo, these indices that you link (supports kallisto versions >= 0.50.1), so use this version if you want to use the index.
If that does not work then I would open an issue at their GitHub.


text:
MMSeqs taxonomy running for over a day
I've been trying to run mmseqs2 on a few metagenomic assemblies and despite my best efforts in reading the wiki and playing with parameters, the process is taking over a day.
In their paper they claim to be 2 - 18x faster than CAT. I ran CAT on the same datasets and they took at most a few hours. So I can only assume that the long run time for mmseqs is down to user error (i.e. me).
I am currently running the commands (general format):
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24

Could someone please let me know where I'm going wrong or what parameters I can tweak to get this running in the expected time frame?

Have you contacted the developers about this? Do any of these issues match the problem you're having?

Hi there! I actually got in running in ~4 hours but then the results are split into 24 semi-corrupted files (I used 24 threads). I have since contacted the developers and will post a solution once I get a response and figure it out :)

After much googling about I realised that I was missing a step! Here are the commands I used to get it running:
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24 --split-memory-limit 500G
mmseqs createtsv metagenome1_DB metagenome1_mmseqs_taxonomy.txt metagenome1_mmseqs_taxonomy_final.txt

The split-memory-limit flag sped up the process and the final line put all the split files back together

Thanks @Rainman


text:
How to identify phages from Phaster in PhageDB
I am running searches for genes using Phaster {PHASTER (PHAge Search Tool Enhanced Release) is a significant upgrade to the popular PHAST web server for the rapid identification and annotation of prophage sequences within bacterial genomes and plasmids.   } and would like to look up the phages that have them on PhageDB The Actinobacteriophage Database to learn more about their distribution. The blast search returns a lot of information about the phages and I am unsure what part of the sequence it gives for phages would be helpful.

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.

What is Phaster? Can you please provide a brief description about the program, and ideally a link that we can visit to get more information? Also... this question doesn't really have a question in it. Could you please clarify what more precisely what you are struggling with?

I am not sure is a good idea to look for result on PHAGEDB , its full name list https://phagesdb.org/data/?set=all&type=simple ,  is probably smaller the the hits out of PHASTER

unless I am getting it wrong and you are asking to use https://phagesdb.org/blast/ local blast to search for your sequence inside PHAGEDB ?

mmmmmh  the phages on PHAGEDB that are not deposited at NCBI are the only one you should look/screen for , but still I am not getting your question PHASTER is used to track leftover presence of phage in plasmid and genomes , PHAGEDB is a repo (kind of non official for actionobacteria phages) or am I wrong ?


text:
How can I convert a BED file to GTF/GFF with gene_ids?
Given a .bed (BED12), how can I convert it to GTF/GFF formats with gene_id attributes? What is the fastest way or available tools to do it?
For example, given an input like this:
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

How can I convert it to:
chr27 bed2gtf gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 bed2gtf transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 bed2gtf exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

Here I provide an ordered list of options (note that I am the author of bed2gtf and bed2gff):
bed2gtf
A high-performance BED-to-GTF converter written in Rust from https://github.com/alejandrogzi/bed2gtf.
Usage: bed2gtf[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gtf)

The isoforms file specification:

a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

Converts

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 3.25 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 1.99 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 1.20 seconds.
Gallus galus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.36 seconds.

bed2gff
A Rust BED-to-GFF3 translator that runs in parallel from https://github.com/alejandrogzi/bed2gff.
Usage: bed2gff[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gff)

The isoforms file specification:

a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

Convert

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 4.16 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 2.15 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 1.30 seconds.
Gallus gallus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.51 seconds.

bedToGenePred + genePredToGtf + refTable
UCSC offers a fast way to convert BED into GTF files through KentUtils or specific binaries using:
bedToGenePred in.bed /dev/stdout | genePredToGtf file /dev/stdin out.gtf

You can install these tools with bioconda, or download them here. The gene_id is only achieved when using refTables (a format specified in UCSC's web browser), you can see a more elaborate answer here Obtaining Ucsc Tables Via Ftp And Converting Them To Proper Gff3 Via Genepredtogtf?.
Other options
Other scripts/tools That DO NOT produce a complete GTF file (lacking gene_id attributes) are:

gtf2bed

gtf2bed < foo.gtf | sort-bed - > foo.bed
awk '{print $1"\t"$7"\t"$8"\t"($2+1)"\t"$3"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' foo.bed > foo_from_gtf2bed.gtf

-kscript
from https://github.com/holgerbrandl/kscript:
kscript https://git.io/vbJ4B my.bed > my.gtf

pfurio/bed2gtf

from https://github.com/pfurio/bed2gtf:
python bed2gtf [options] <mandatory>

AGAT

AGAT
Considering only the options that produce gene_ids attributes, bed2gtf and bed2gff are faster by ~3-4 seconds than UCSC's C binaries. More detailed instructions of this tools are explained in the sources linked.

Thank you. These are great resources.


text:
single_to_multi_fast5 does not collect all the single files if the input folder contains mixed types of fast5 files
I have a dataset that contains thousands of mixed multiple and single fast5 files in a non-homogenous folder structure.
The reads and fast5 files are coming from different labs and some labs gave single fast5 files and some multi fast5 files, without clear mappings. My final goal is to convert all the fast5 files to multi fast5 files.
There are commands in ont_fast5_api which is an interface to HDF5 files of the Oxford Nanopore.
Using this API we can convert single fast5 files and multi fast5 files to each other.
this is the link to this API github page:
https://github.com/nanoporetech/ont_fast5_api
So I made a test mixed folder to check if I can use the API without missing reads as a result of conversion. The test folder structure is similar to my original folder structure:

mixed

Subfolder_A

single_fast5 x 4000

Subfolder_B

multiple_fast5 (contains 4000 single reads fast5)

single_fast5 x 2
multiple_fast5_1 (contains 4000 single reads fast5)
multiple_fast5_2 (contains 4000 single reads fast5)
multiple_fast5_3 (contains 100 single reads fast5)

Here is a screenshot of the structure with a few files from each folder:

and I created two other folders outside of the mixed folder as follows:

single

multi

Here my question is about single_to_multi_fast5 and multi_to_single_fast5 commands.
test_path=mixed
multi_path=multi
single_path=single
multi_to_single_fast5 -i $test_path/ -s $single_path/ --recursive

in the $single_path, I expect to see, all of the reads in multi files converted to single files (12100 reads in 12100 single fast5 files), but I see 8100 single files.

So I have 4000 missing reads.
Now I convert single read fast5 file to multi:
single_to_multi_fast5 -i $test_path/ -s $multi_path/ --filename_base $output_name --batch_size 1000 --recursive

in the $multi_path folder, I expect to see that generated multi files contain all my single-read files (4000 + 2), but I see 3006 in the mapping summary.

So again I missed 996 single reads in my multi files.
Since I made a small subsample of data, I am sure that the files are not overwritten.
Because my final goal is converting everything to multi I was going to do the following pipeline:
orig_path=original_folder
intermediate_path=intermediate
final_path=final_multi
# converts/collects all multi files in the original folder
multi_to_single_fast5 -i $orig_path/ -s $intermediate_path/ --recursive
single_to_multi_fast5 -i $intermediate_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

# converts/collects all single files in the original folder
single_to_multi_fast5 -i $orig_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

However because it seems I get missing reads during the process, I cannot use this pipeline.
At the same time, my original dataset is super huge, I cannot check if individual files are single or multi. It would take ages (I guess).
Is there a solution for solving this problem besides checking every file? Maybe my pipeline is not quite suitable for what I want to do.

Can you please show the first few filenames in each folder? e.g. by running the following command: find mixed -type 'd' | while read dname; do echo "** ${dname} **"; tree -L 1 ${dname} | head; done

Thank you for your reply! I included the first few filenames in each folder @gringer

The image seems like it's annotated incorrectly. A file with a channel and read number should be a single fast5 file, not a multi fast5 file.

Yes, the 2nd label (multi fast5 files with 4000 reads) and the 3rd label (single fast5 files) should be swapped. @gringer

Okay, that's good to know; thanks for the clarification.

It looks like you're creating too much additional work for yourself. What I'm trying to lead to with my comments is that the file names should help in working out if a file represents single reads or multiple reads, in which case the single reads can be shifted into a separate folder.
Based on the file names that I would expect, a search for files containing '_read_[0-9]' should be sufficient to fish out the single fast5 files:
find . -name '*_read_[0-9]*' > single_reads.txt
mkdir -p single_reads
cat single_reads.txt | while read fname; do mv -i ${fname} single_reads; done

[That single_reads.txt file could be split into multiple files if directory entries get too large.]
After which you should have a single folder that contains thousands of single-read fast5 files, and other folders that contain multi-fast5 files, removing the problem of needing to deal with mixed file types.

I am gonna try this, thank you! @gringer


text:
BioPython bootstrap is not reliable?
Here i will show you a minimal working example of code and as you can see the support values for the tree is always 100.
I am using synthetic sequences of 100bp for 6 elements. The sequences have been generated at random choosing from ATCG for each position with equal probabilities.
The input fasta file is
>S1
CTCAGGCACTAGGAGCTTCTCCAGGGCAAAGTTGTCTACGAATATGCCGACTCAGAAGGTTATCAATACGGTTTACTTATCTGCACGCCATTTTCCTATG
>S2
TAACAAATGTTTCTCGCTCGAACCCCGTGGTGCGAGGGGTCACGATAAAAGGTCCCTCCTTCCACGGCATAATTGCTCCCTTTCTTTTCCGTGGGCAGGA
>S3
TTCCTAAAGGACGACTGGAAGCCGGGTACACGTCAACGGAGCTTCTAGCCCGAGGTCTATAGACCGGTTAATGAACTAGACTTAATGTGGAGCTCGTAGA
>S4
ACTTCGTGCCCAACCACTCTATGAGCAGAGTGCTCGAATGAACCTCAAAAGGATATCCGCATTTACTCTATATAACAAACGGCCGTCCCCCCATCTGTCA
>S5
TGACACTGGGTCATTTTACCCGCACTGATCGCTGGGCAGGTCGGCAATTTCGTCAGAATGCCGTGGCCGCACTGAAAATATCATTACCGGACTGAGTATG
>S6
ATAAGCCGAGGGGTAGCCCTCTATTTCGCACCGTATAGACGAAGTGATAAACTTTCTAACCACGTGTCGCCTATCTCTACCTAGCCACATTTGAGGTGCG

The following python code makes the multiple alignment using MUSCLE algorithm
### alignment
from Bio.Align.Applications import MuscleCommandline
from Bio import AlignIO, SeqIO
import os

muscle_exe = "C:/Scripts/muscle3.exe"
path = os.getcwd()
with open('tmp_in.fa', 'w') as fp:
    fp.write(fasta)

## replacing \\ with /
inp_fp = path.replace('\\', '/') + '/tmp_in.fa'
out_fp = path.replace('\\', '/') + '/tmp_out.afa'

muscle_cline = MuscleCommandline(muscle_exe, input=inp_fp, out=out_fp)

The following code block takes the previous alignment out_fp and creates a tree using UPGMA algorithm.
The number of bootstrap trees is set to 50 and the consensus is majority based.
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo.Consensus import *

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

The final consensus tree looks like this
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.350925925925926, name='S5')
        Clade(branch_length=0.05833333333333341, confidence=100.0)
            Clade(branch_length=0.05879629629629617, confidence=100.0)
                Clade(branch_length=0.26481481481481484, name='S3')
                Clade(branch_length=0.020370370370370372, confidence=100.0)
                    Clade(branch_length=0.24444444444444446, name='S1')
                    Clade(branch_length=0.24444444444444446, name='S4')
            Clade(branch_length=0.031018518518518404, confidence=100.0)
                Clade(branch_length=0.2925925925925926, name='S6')
                Clade(branch_length=0.2925925925925926, name='S2')

My question is: why does the tree have all support values equal to 100? I am sure there must be an error either on my side or in the package module itself.

I think this is a bug.
It seems to work if you do this, creating an equivalent Alignment object instead of a MultipleSeqAlignment to give the bootstrap step:
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)

bootstrap_consensus calls bootstrap_trees, which makes however many randomly shuffled alignments you asked for.  But if the input is a MultipleSeqAlignment object like you're using it builds a tree over and over with that same input to the function as the tree function input, so it always gets the same tree over and over.  The other branch of bootstrap_trees, which runs if you have an Aignment object, builds trees based on the randomly-selected columns of the input, like you'd expect.  Pretty sure that's just a variable naming screw-up one of us should report to Biopython.
Self-contained version of your tree part:
#!/usr/bin/env python
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio import AlignIO
from Bio.Phylo.Consensus import *

out_fp = "tmp_out.fa"

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

# Give bootstrap_consensus an Alignment object instead of MultipleSeqAlignment
# object
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)

Example outputs for the two different input alignment objects:
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.012239583333333349, confidence=100.0)
            Clade(branch_length=0.31666666666666665, name='S3')
            Clade(branch_length=0.033333333333333284, confidence=100.0)
                Clade(branch_length=0.2833333333333333, name='S2')
                Clade(branch_length=0.2833333333333333, name='S4')
        Clade(branch_length=0.033072916666666716, confidence=100.0)
            Clade(branch_length=0.29583333333333334, name='S5')
            Clade(branch_length=0.03750000000000007, confidence=100.0)
                Clade(branch_length=0.2583333333333333, name='S1')
                Clade(branch_length=0.2583333333333333, name='S6')
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.2583333333333333, name='S2')
        Clade(branch_length=0.03405145202020203, confidence=66.0)
            Clade(branch_length=0.275, name='S5')
            Clade(branch_length=0.049892241379310344, confidence=58.0)
                Clade(branch_length=0.21250000000000002, name='S1')
                Clade(branch_length=0.21250000000000002, name='S6')
        Clade(branch_length=0.039416666666666676, confidence=50.0)
            Clade(branch_length=0.2583333333333333, name='S4')
            Clade(branch_length=0.321875, name='S3')

Thank you very much i think this was the issue. I will report it to their github if you haven't already raised an issue.

I haven't yet, so, go for it!  Should probably be a simple fix for them.

Initially, 100 looks like a good value.  However, if you lower to bootstrap value, let's say maybe to 50 as you did, you will probably get an erroneous support value.  Bootstrapping (you probably already know this) is how many times a branch in a phylogenetic tree is observed.  A good bootstrapping value would be anywhere between 80 and 100.  50 is not a good bootstrapping value.  This means the branches are only being observed 50 times as opposed to 85, 90, etc...   If I were you, try changing the bootstrapping value of 50 to at least 80.  See what that gets you.  You may be satisfied with it, but you may want to increase it depending on your output.  By doing this, you will see a more realistic and phylogenetically accurate output of your tree.  Most bioinformaticians use 100 as a bootstrap value.  The only issue, I wouldn't really call it an issue, but it will take longer for an output to occur due to computing time and computing memory available.  BUT you will still get an output, you just may need to wait a bit.   My bioinformatics background is rooted in phylogenetic and horizontal gene transfer research involving phages and bacteria, so I have seen this issue a lot and have experienced myself a lot!  I hope this information provides some insight to your project and bioinformatics journey! :)

I don't think this is the case. This is a Minimal working example to show the problem. The real data is hundreds of species with longer sequences. When I run multiple times manually the clades in the tree change, but when I use biopython consensus tree it does not happen, so it must definitely be a bug or an error in my code.

When you mean you ran the samples "multiple times manually" as opposed to "when I use the biopython consensus tree"?  This sentence is confusing and misleading potentially.  "Consensus trees are used to summarize a set of trees defined on the same set of taxa. A consensus algorithm takes the trees as input, so that the method of producing the input trees is not part of the consensus algorithm" as cited on:  https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2909780/#:~:text=Consensus%20trees%20are%20used%20to,part%20of%20the%20consensus%20algorithm.

A consensus tree uses the consensus algorithm to generate a tree.   That makes me wonder if the algorithm you are using is appropriate for the amount of samples you have and the variation of species from different axons that you are using.  What taxons are you using?  Are you looking at family, class, order, genus, etc...?  Family would defiantly be a broader range of values than genus.  If you can maybe narrow down to a more exclusive taxon like order or genus, maybe your values may be more consistent?

In the example i gave (that is completely sufficient to understand the problem), there is no mention of taxa. There is no need to invoke more complex processes when it is evident that when you run a consensus tree either on a very random sequence set (MWE) or in a very big dataset (my actual case) you do not expect to see 100% on all the nodes. It's like saying that you made the best possible alignment, also all the assumptions of your models are met and there is no measure error in the sequences and so on...