Hey everyone,

I’m currently working on a Python script using the Bioformats library to process .lif files. My goal is to extract everything contained in these files (images and .xml metadata), essentially replicating what the Leica software does when exporting data.

So far, I’ve managed to extract all the images, and at first glance, they look identical. However, when comparing pixel by pixel, they are actually different. I suspect this is because the Leica software applies a LUT (Look-Up Table) transformation to the images, and I haven't accounted for that in my extraction.

Another issue I’m facing is the .xml metadata file. The one I generate is completely different from what Leica produces, and I can’t figure out what I’m missing.

Has anyone encountered a similar issue? Does Bioformats handle LUTs differently, or should I be using another library? Any suggestions on how to properly extract the correct images and metadata?

I’d really appreciate any insights! Thanks in advance.

Slurm is everywhere in academia, but what about biotech and pharma? A lot of companies lean on cloud-based orchestration—Kubernetes, AWS Batch, Nextflow Tower (I still think they're too technical for end users)—but are there cases where Slurm still makes sense? Hybrid setups? Cost-sensitive workloads?

If you work (or have worked) in private-sector bioinformatics, did Slurm factor into your workflow, or was it all cloud-native? Curious what’s actually happening vs. what people assume.

I’m building an open-source cluster compute package that’s like a 100x simpler version of Slurm, and I’m trying to figure out if I should just focus on academia or if there are real use cases in private-sector bioinformatics too. Any and all info on this topic is appreciated.

Hi everyone,

I'm trying to cluster the proteomes of 477 P. aeruginosa into orthologs and having some difficulty with Orthofinder. Initially running it on all 477 took far to long to compute on our cluster, so I selected a core of 15 which have the phenotypic traits I am interested in. I then added in the rest of the species with the --assign option.

Out of 2939270 genes, this has resulted in 11174 not being assigned to orthogroups (0.38%). After refining this to HOGs, an extra 5922 are then not placed into any HOG at the N0 level. Whilst this is a small fraction of my dataset, I'm unsure why this is even happening at all. I've checked the Orthogroups_UnassignedGenes file, but that only contains 183 genes and all of them are assigned to orthogroups anyway, just orthogroups with a size of 1. These genes aren't limited to any particular bacteria, with 389/477 having at least one gene not in an orthogroup. The number unassigned genes ranges from 1 - 425.

Does anyone have any insight on why this could be occurring? I've opened an issue on the github page but the developers don't seem to be super active with their latest response being over 3 weeks ago. I'm not even sure on the best thing to do next to troubleshoot!

Thanks in advance

Hi everyone, cross posting from r/MachineLearning in the event this is a better venue. I’ve been working a custom dataset that I’ve curated and I’m about worn out with it:

I’m working on a hobby project in which I’ve collected complete proteome sequences for flu isolates collected around the world from about the year 2000 to the present. As you can imagine, this real world data is plagued with recency bias in the number of isolates recorded, and their are many small minor classes in the data as well (single instance clades for example).

For context, there are many examples in the literature of modeling viral sequences with a variety of techniques, but these studies typically only focus on one or two of the 10 major protein products of the virus (Hemagglutinin (HA) and Neuraminidase (NA)). My goal was to model all 10 of these proteins at once in order to uncover both intra- and inter- protein interactions and relationships, and clearly identify the amino acid residues that are most important for making predictions.

I’ve extracted ESM embeddings for all of these protein sequences with the 150M param model and I initially trained a multi-layered perceptron classifier to do multi-task learning and classification of the isolates (sequence -> predict host, subtype, clade). That MLP achieved about 96% accuracy.

Encouraged by this, I then attempted to build predictive sequence models using transformer blocks, VAEs, and GANs. I also attempted a fine-tuning of TAPE with this data, all of which failed to converge.

My gut tells me that I should think more about feature engineering before attempting to train additional models, but I’d love to hear the communities thoughts on this project and any helpful insights that you might have.

Im currently working on doing molecular docking and screening for drug discovery. Currently I have screened ~2000 drugs and I have compiled all their ADMET properties. I was wondering how I should rank the bioavailability for each drugs based on their ADMET subcategories scores. Is there like a ranking of what subcategories are more important compared to others?

Here are all the properties that were screened for.

Caco_2_c Caco_2 HIA MDCK F50 F30 F20 BBB OATP1B1_inhibitor OATP1B3_inhibitor OATP2B1_inhibitor OCT1_inhibitor OCT2_inhibitor BCRP_inhibitor BSEP_inhibitor MATE1_inhibitor Pgp_inhibitor Pgp_substrate PPB VDss CYP1A2_inhibitor CYP3A4_inhibitor CYP2B6_inhibitor CYP2C9_inhibitor CYP2C19_inhibitor CYP2D6_inhibitor CYP1A2_substrate CYP3A4_substrate CYP2B6_substrate CYP2C9_substrate CYP2C19_substrate CYP2D6_substrate HLM RLM UGT_substrate CLp_c CLr T50 MRT Neurotoxicity DILI hERG_1uM hERG_10uM hERG_30uM hERG_1_10uM hERG_10_30uM Respiratory_toxicity Nephrotoxicity Skin_sensitisation ADT Ames Mouse_carcinogenicity_c Mouse_carcinogenicity Rat_carcinogenicity_c Rat_carcinogenicity Rodents_carcinogenicity Micronucleus Reproductive_toxicity Mitochondrial_toxicity Hemolytic_toxicity Repeated_dose_toxicity AOT_c AOT FDAMDD_c FDAMDD AR ER AR_LBD ER_LBD Aromatase AhR ARE ATAD5 HSE p53 PPAR MMP TR GR

Hi all,

I have 18 genes of interest in a reference strain of bacteria which are all next to one another. I would like to see if they are all conserved in my other isolates (n=11) and in the same order.

They are not at the same coordinates as the assemblies are not rotated to dnaA and do not have the same locus ID's because PGAP doesn't seem to keep them consistent between genomes.

My aim is to draw a gene arrow plot in gggenes to visulise the suspected rearrangements. Is there a quick way to pull the genes out of a multi-fasta or similar file and make this all work?

EDIT: example of the figure i'm trying to achieve

https://preview.redd.it/w6xw66oxrage1.png?width=576&format=png&auto=webp&s=db79af4b17ae15db615d81dc9443d6427ea906ec

Hi, I am trying to do Transcriptome analysis with the RNAseq data (I don't have bioinformatics background, I am learning and trying to perform the analysis with my lab generated Data).

I have tried to align data using tools - HISAT2, STAR, Bowtie and Kallisto (also tried different different reference genome but the result is similar). The alignment score of HIsat2 and star is awful (less than 10%), Bowtie (less than 40%). Kallisto is 40 to 42% for different samples. I don't understand if my data has some issue or I am making some mistake. and if kallisto is giving 40% score, can I go ahead with the work based on that? Can anyone help please.

I’m considering using an unsupervised clustering method such as kmeans to group a cohort of patients by a small number of clinical biomarkers. I know that biologically, there would be 3 or 4 interesting clusters to look at, based on possible combinations of these biomarkers. But any statistic I use for determining starting number of clusters (silhouette/wss) suggests 2 clusters as optimal.

I guess my question is whether it would be ok to use a starting number of clusters based on a priori knowledge rather than this optimal number.

I'm working on a pathogen that is able to infect a lot of hosts and is also present around the world. I have sequences from samples around the world and also from 53 different hosts. Building the tree is easy, but is there a software/package that lets me easily add extra info like the host and location so I can spot patterns or will I have to do that manually?

Most of the questions I've seen about batch correction were about integrating samples from different datasets (already annotated). I wanted to get some opinions about whether the best practice would be to annotate each sample individually then integrate + batch correct or integrate and then annotate the object together after batch correction. Maybe it doesn't matter with unannotated data (newly generated datasets from the lab or collaborators)?

I currently do the latter but maybe there's some disagreements with this? Just wanted to double check i'm not botching the whole scRNA pipeline with what I'm doing... lol

> head(colnames(seurat_object)) [1] "JC_Dx_AAACCTGGTCGCTTCT-1" "JC_Dx_AAACGGGAGGTAGCCA-1" "JC_Dx_AAACGGGGTAAAGGAG-1" [4] "JC_Dx_AAAGATGCACTGTCGG-1" "JC_Dx_AAAGATGGTCTGCCAG-1" "JC_Dx_AAAGATGTCCTGTAGA-1" > head(combined_TCR$JC_Dx$barcode) [1] "JC_Dx_AAACCTGAGTACGACG-1" "JC_Dx_AAACCTGCAACACGCC-1" "JC_Dx_AAACCTGCAGGCGATA-1" [4] "JC_Dx_AAACCTGCATGAGCGA-1" "JC_Dx_AAACGGGAGCGTTTAC-1" "JC_Dx_AAACGGGAGGGCACTA-1"

Hi I converted my adata object to a seurat object using createSeuratObject and have the colnames of my seurat object that look with the sample_id_barcode-1 which is the same format my combined.TCR object barcode looks with the sample_id as the prefix and when I use the combineExpression(combined.TCR,seurat_object,group_by="sample") gives me an error:
Warning message: In combineExpression(combined_TCR, seurat_object, cloneCall = "aa", : < 1% of barcodes match: Ensure the barcodes in the single-cell object match the barcodes in the combined immune receptor output from scRepertoire. I even tried the method in the faq and nothing seems to help. Any help would be greatly appreciated

I've found the posts about samtools and the other applications that can accomplish this, but is there anywhere I can get this done without all of those extra steps? I'm willing to pay at this point.. I have a CRAM and crai file from Probably Genetic/Variantyx and I'd like the VCF. I've tried gatk and samtools about a million times have no idea what I'm doing at all.. lol

How do I split my .bam into 1000nt bins using Mac terminal. Thanks.

I am working on a simple fastq -> mutant enrichment score pipeline, but wonder if I'm not thiking to simplistic. This is the idea...

Setup:

• I have an UNSORTED and SORTED sample, 2 fastqs each.. R1 and R2. Readlenght is 150bp.
• The sequence of interest is a 192bp long sequence.
• R1 has a primer1 indicating the start of sequence of interest
• R2 has a primer2 indicating the start of sequence of interest

My approach

1. Trim raw data using the primers, keeping only the region of interest
2. Merge R1 and R2, creating the complete region of interest (discarding all resulting reads not being 192bp and filtering on quality 30). Little of over 80% of reads remain here btw.
3. (Use seqtk to) translate DNA sequence to protein sequence (first fastq to fasta, then fasta to protein)
4. Calculate frequency of protein mutants/variants (nr of variants divided by total amount) for each sample
5. Calculate enrichment using ratios from 4) (freq-SORT/freq-UNSORTED)?
6. log2 transform the results from 5)

End result:

Data table with amino acids sequence of interest as cols, amino-acid changes as rows and log2(enrichmentratios) as values which will then be plotted in the form of a heatmap based on enrichment ratios...

Because we are looking at a fixed sized sequence which is entirely within the PE reads no mapping is necessary.

I have been looking into various options for DMS (enrich2, dms_tools2, mutscan) but if the above is correct then diving into those tools feels a bit much...

I feel like I'm looking at iit too easy though, what am I missing?

*EDIT

We have been able to compare the results from this with earlier generated data and even though the exact enrichment values matter, the trend (enrichment) is just about perfectly overlapping... So still looking into what we might be missing but at least the approach corresponds to what was done before

Hi everyone, I've read an article where they built a database includes about 10k molecules and calculate the TCs distribution of all (based on 1024bit ECFP4 ). It doesn't develop their own way to calculate it but cites a method from a paper published in 2000 and the SVL code used is not avalible anymore. So I googled it and only find this one but this program is also obsolete.

So I wonder which program/software might gives this function? Maybe they self-built a complex program and executed this calculation completely in RDkit?

Hi there. I am working with metagenomics data and I am using the MBECS package to perform batch correction on the data. I have 2 batches (both done on the same MiSeq sequencer), one with 6 samples and one with 74 samples (both with 50% cases and controls aprox.).

I have used Principal Least Squares Discriminant Analysis (PSLDA) as method for the batch correction.

After applying the batch effect correction, I see a reduction on the batch effect according with the follwing Principal Variance Component Analysis (PCVA). Raw clr-norm data is represented on the right and PSLDA batch-corrected data in on the left.

Principal Variance Component Analysis (PCVA). Left: Uncorrected data. Right: Batch-corrected data.

Nevertheless, despite the seq_batch (sequencing batch) explained variance goes down to 0%, the interaction between batch and group increases by ~3X.

Can someone explain why does this happens? Shouldn't it be reduced since batch effect is corrected?

Also looking at the PCA, seems that the batches are now more clearly separated after batch correction, but from the other side, silhouette coefficient shows less difference between bathes.

Principal Component Analysis (PCA). Top: Uncorrected data. Bottom: Batch-corrected data.

Silhouette Coefficient. Left: Uncorrected data, Right: Batch-corrected data.

Can anyone throw some light on this? Do you think is worth it to apply batch correction?

Thank you very much in advance.

I am currently analyzing RNA-Seq data from human samples. The sequencing was done by Novogene using an lncRNA library preparation (not polyA-enriched).

I aligned the raw reads to the latest human reference genome (Ensembl) using HISAT2, achieving >90% mapping rates for all samples. However, when quantifying mapped reads using featureCounts, I observe that the assigned reads are much lower—ranging from 30% to 55%.

I am trying to understand whether this is a technical issue or expected due to the higher sequencing depth (~12 Gb per sample) and the lack of polyA enrichment.

Here this the code I used:

featureCounts -a "$GTF_FILE" -o "$output_file" -p -T 16 $bam_files -g gene_id --countReadPairs -s 2

Any input on this will be greatly appreciated!

Hello. I am hoping this is the correct place to post and apologies in advance for maybe using the wrong terminology. I am currently a masters student studying mathematics, and for my dissertation I am looking at applying graph invariants to biological networks. My plan is to start on smaller networks so that I can do some calculations by hand but I am having a difficult time finding appropriate networks or being able to understand what I am being shown. I am using STRING database and have somewhat figured out how to tailor it to what I am looking for but my question is, say in the image I have uploaded, STRING is telling me that there are 6 edges, which I can see obviously. However, I do not understand what the different colours represent and if that is relevant, if I am looking at networks in a mathematical sense rather than a biological. If they are relevant, how is the best way to go about understanding this more? Again, apologies if my question isn't clear, this is all very new to me. Thank you for any advice/help you can offer.

https://preview.redd.it/4koa3ijw51ge1.png?width=3250&format=png&auto=webp&s=b417f3f41e28f5b41ac6bed357c509b5d551c2bc

What is a normal expected alignment rate for cfDNA onto a reference genome? My data is cfDNA mapping onto a mouse genome (mm39), but literally any number with a citation will do. I'm having a very difficult time finding a paper that reports an alignment rate for cfDNA onto a reference genome, and I just want to know what is an expected range. Thanks !

I'm using BWA MEM as an aligner, but it could be another as well.

My background is in computer science and I recently started developing software and pipelines for bioinformatics... I'm from Brazil and bioinformatics software/pipeline development is not common field here. Because of this, I'm also considering submitting to an international journal (open access) to get more visibility. Any suggestions?

Hello, when we have a dataset of Single cell RNA-seq of a given cancer type in different stages of development, do we utilize a supervised analysis or unsupervised approach?

I am working from a Drosophila dm3 gtf file trying to infer different percentage compositions of genomic features of interest (UTRs, CDS, introns, etc.) Since there is no "intron" feature explicitly found in the file I decided to obtain them by:

1. bedtools merge on file only containing "transcripts"

2. bedtools merge on file containing the remaining features (CDS, exons, UTRs, start, and stop codons)

3. bedtools subtract using - a "transcripts" file and -b "remaining_features" file

4. Use awk '{total += $3 - $2} END {print total}' intron_file.txt to calculate total intron bp

5. Total intron bp / Total Drosophila dm3 genome bp where total genome bp was obtained from (https://genome.ucsc.edu/cgi-bin/hgTracks?db=dm3&chromInfoPage=)

The value I get is usually >42% compared to the 30% mentioned in literature (Table 2 from Alexander, R. P., Fang, G., Rozowsky, J., Snyder, M., & Gerstein, M. B. (2010). Annotating non-coding regions of the genome. Nature Reviews Genetics, 11(8), 559-571. )

What could I be doing wrong? Things I should look out for? Thank you for the help!

I'm on a local system and i recently learned many new bioinformatic techniques and discovered new pipelines which I would like to try and test out myself on some data. However,I'm on a local system and not on the cluster and I am looking for a platform to try out the various codes and analyses based on open source data from previous publications (and essentially retrace the results). This coming month I'm willing to try Mostly some ATAC seq pipelines, and snp calling with some RNAseq if I have time. I'm a novice in this matter. Please feel free to give as many inputs as you want.

I had a question on determining when to use each of these sequencing methods. Asked my prof about this but he wasn’t very clear on it :/ Also, when conductinf paired end readings with Sanger, are the paired end reads done by pcr, then another subsequent sequence via sanger? Or are they done in one round?

Thank you!

I am in the process of down-sampling 10x multiome data (paired scRNA and scATAC) due to differences in depth per cell of final libraries and I am trying to determine which FASTQ files to down-sample for the ATAC portion. It looks as though the samples contain dual indexing and as such, each sample has an R1, R2, I1, and an R3 fastq file. 
According to the 10x website here the I1 and R2 reads contain indexing information. Is it correct to down-sample the R1 and R3 fastq files or do the indexing files also need to be downsampled?

Currently doing this with Seqtk specifying a consistent random seed. GEX went smooth but really not sure how to handle the ATAC portion.

Has anyone ever tried using the downsampleReads function from DropletUtils R package to achieve this in a less cumbersome way? I know it will work fine for the GEX portion, but not sure how it will handle the ATAC.

I am analyzing a pbmc/tumor experiment

In the general populations(looking at the oxygen groups) the CD14 dot is purple(high average expression) in normoxia, but specifically in macrophage population it is gray(low average expression).

So my question is why is this? Because when we look to the feature plot, it looks like CD14 is mostly expressed only in macrophages.

This is my code for the Oxygen population (so all celltypes):

Idents(OC) <- "Oxygen" seurat_subset <- subset(x = OC, idents = c("Physoxia"), invert = TRUE)

DotPlot(seurat_subset, features = c("CD14"))

This is my code for the Macrophage Oxygen population:

subset_macrophage <- subset(OC, idents = "Macrophages") > subset(Oxygen %in% c("Hypoxia", "Normoxia"))

DotPlot(subset_macrophage, features = c("CD14"), split.by = "Oxygen")

Am i making a mistake by saying split by oxygen here instead of group by?

There an ongoing fuss about deepseek . Has anyone tried it to try provide code for a complex bioinformatics run and see how it performs?

I’m wondering—are there people working in bioinformatics who use Rust? Most tools seem to be written in Python, C, or R, but Rust has great performance and memory safety, which feels like it could be useful.

If you’re in bioinformatics, have you tried Rust for anything?

Hi all,

I have a relatively large set of sequence logos for a protein binding site. I am interested in comparing these (ideally pairwise). Trouble is, I haven't been able to find much as far as metrics to compare sequence logos. In my imagination, I would like something to the effect of a multi-sequence alignment of the logos, from which I then have a distance metric for downstream analyses. The biggest concern I have is the compute time that could be required to make all of the comparisons. Worst case scenario, I will just generate an alignment with the ambiguous strings. Alternatively, I will fix the logo size and could try to come up with a method to determine edit distance between these strings.

One final (probably important detail) is that I am working with nucleotide data and looking at logos between 8-16 base pairs.

Any help is definitely appreciated!

https://preview.redd.it/xupe0hnaevfe1.png?width=1385&format=png&auto=webp&s=665c31834ce3a2ff05ab4f04ef3df63d3f1c8b7b

This is a figure from analysis of scMultiome dataset (https://doi.org/10.1126/sciadv.adg3754) where the authors have shown the concordance of RNA and ATAC clusters. I am also analyzing our own dataset and number of clusters in ATAC assay is less than RNA, which is expected owing to sparse nature of ATAC count matrix. I feel like the figure in panel C is a good way to represent the concordance of the clusters forming in the two assays. Does anyone know how to generate these?