text:
Warning message related to tensorflow package when trying to run deepbgc bioinformatics tool
I have come across some related issues pertaining to this tensorflow warning but none related specifically to deepbgc
I have created a conda environment and have installed onto it the bioinformatics analysis tool deepbgc. Now the program runs but it does output some warning messages related to the tensorflow package. When I run the "deepbgc info" command through terminal, I get the following output messages related to the tensor flow package:
Using TensorFlow backend.
WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:74: The name tf.get_default_graph is deprecated. Please use tf.compat.v1.get_default_graph instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:517: The name tf.placeholder is deprecated. Please use tf.compat.v1.placeholder instead.

WARNING 03/10 12:55:05   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:4138: The name tf.random_uniform is deprecated. Please use tf.random.uniform instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:133: The name tf.placeholder_with_default is deprecated. Please use tf.compat.v1.placeholder_with_default instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:3445: calling dropout (from tensorflow.python.ops.nn_ops) with keep_prob is deprecated and will be removed in a future version.
Instructions for updating:
Please use `rate` instead of `keep_prob`. Rate should be set to `rate = 1 - keep_prob`.
WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:174: The name tf.get_default_session is deprecated. Please use tf.compat.v1.get_default_session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:181: The name tf.ConfigProto is deprecated. Please use tf.compat.v1.ConfigProto instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:186: The name tf.Session is deprecated. Please use tf.compat.v1.Session instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:190: The name tf.global_variables is deprecated. Please use tf.compat.v1.global_variables instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:199: The name tf.is_variable_initialized is deprecated. Please use tf.compat.v1.is_variable_initialized instead.

WARNING 03/10 12:55:06   From /EFS/tools/miniconda/envs/deepbgc2/lib/python3.7/site-packages/keras/backend/tensorflow_backend.py:206: The name tf.variables_initializer is deprecated. Please use tf.compat.v1.variables_initializer instead.

Could this be related to the version of tensorflow? Currently, I have version 1.15.4 installed onto the conda environment. Are these error messages significant and could they impact the data running through deepbgc? Or can they mostly be ignored? Any guidance appreciated, thanks!

That should be easy to sort out and its good to get it sorted. I'll reply tomorrow.

@M_ great, I appreciate it. I was thinking maybe this was an issue with my tensor flow version so I tried updating it but the problem with that is it makes other packages like numpy incompatible with the version of deepbgc I'm running on. Could an update to deepbgc do the trick and then subsequently, an update to tensorflow?

The 'blunderbus' way to do this is as follows (not recommended in this instance :-) ),

    conda update --all
    conda update -n base conda
    conda update -n myenv --all

The refined way to do this is

    conda update -n myenv tensorflow # see `Best guess, notes

In theory conda was supposed to sort out the deprecation, because all the versions were supposed to be harmonious. I'm a bit puzzled (see below). What might have happened is tensorflow has been upgraded (or keras) against the version the developer had written it for.
Here myenv is the name you gave that given conda environment.
In your case, I would look at point 2, because the environment is very specific to your project and ensures reproducibility.
Note A purist would specify the exact channel -c that tensorflow is found in, or whatever dependency. This actually can be very important because there can be name conflicts, obviously not for tensorflow because I'm pretty sure Google have the name registered.
The reason I do point 1 is the problem might a dependency of tensorflow and you'd have to figure it out, but again see 'notes'.

Best guess
In this particular example,
conda update -n myenv keras

Might be the ticket, because what I think is going on is its not actually tensorflow but the keras sitting above it.
The alternative (probably better) guess is ....
conda update -n myenv -c bioconda deepbgc

This is probably what I would try first.
Note I'm surprised that conda/mamba has not resolved this issue on installation because thats exactly what is it is supposed to do. In other words its supposed to make a perfect installation of deepbgc precisely as the developer had intended it. There is a risk that updating might not be the solution ... and solution is a new conda/mamba installation (see note 2).
Note 2 Watch for the Python version, you are using 3.7 ... that might be an issue. If the developer had written this for >3.7 that could cause problems. Different Python versions have different methods and if one of the 'new' methods was being used that would be an issue. Ultimately your  package/pipeline is really complex and the more stuff the more chance stuff doesn't fit together. We're currently on 3.10 and 3.9 is really stable.
BTW 3.11 will be a major upgrade of Python3 and will be the 'goto' version.

I looked at deepbgc its built on Python3.7, I don't know why because it was released in 2019. My advice is start again with a new install and follow the instructions exactly as on Github. The order of the build is really important (omitting pip use conda/mamba instead). It doesn't seem to be very conda but thats how the developer did it. In theory the whole pipeline could have been placed on conda. Instructions here https://github.com/Merck/deepbgc
The key issue is that it is built on an old build of Python, with a lot of dependencies which have substancially changed. There is concern about ...
conda create -n deepbgc python=3.7 hmmer prodigal

because the author is not specifying the versions of hmmer or prodigal. These should have been built into the conda upload.

Thanks! I will go ahead and give this a try. Now just want to confirm this, the conda update -n myenv ... command is run outside the environment, correct? In other words, when the environment is deactivated?

Yes outside the given env. In this case I don't think it will make any difference but its you never know.

So when I try running conda update -n my env keras, I am getting the following error: PackageNotInstalledError: Package is not installed in prefix.
prefix: /EFS/tools/miniconda/envs/deepbgc2
package name: keras - should i be providing the full path to the package when i update? I tried to cd into the directory where the package is located and I got the same error message

I've updated the post


text:
Mapping PDB ID + chain ID to UniProt ID
This question was addressed here.
Unfortunately, including bioservices, the solutions pypdb and map_pdb_to_uniprot do not work for now.
Besides, I do not need residue level mapping as one of the suggested solution.
Does the failure of known methods have anything to deal with the recent update of the UniProt database?
Please let me know if you know a workaround.
Thank you!

https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html

Please find the solution here:
https://www.biostars.org/p/9540519/#9540826

source: Stack Exchange

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

sort: votes


text:
How to promote assemblies into genomes in NCBI?
Note: I've never submitted an assembly/genome to NCBI, so excuse if my perspective is flawed.
I'm working with Drosophila subobscura. (spring fruit fly)
I see here https://www.ncbi.nlm.nih.gov/data-hub/genome/?taxon=7241 that there are at least 2 assemblies of better quality than the only available genome (I checked checked N50 and number of scaffolds).

https://www.ncbi.nlm.nih.gov/data-hub/genome/GCA_903684685.1/
https://www.ncbi.nlm.nih.gov/data-hub/genome/GCA_018903505.1/

However, these assemblies weren't promoted to genomes.
Questions

What is required to consider them genomes?

Is it the scaffold to chromosome notation that is missing? Or also the annotations?

How would I promote the better assemblies into genomes?

And also, if I wanted to annotate one of the better assemblies, how would I proceed?

I'd suggest splitting the fourth question off into its own question. Genome annotation is a very different subject than the NCBI controlled vocabulary and content management procedures. Annotation is also much easier to answer, as it's merely a technical question! E.g. NCBI has their own annotation pipeline, or you could use MAKER.

The reason is simple, insect genomes are often heavily fragmented and historically relied on wet-lab experimentation - polytene chromosome hybridisation to correcty assemble the chromosomes. Thus you probe the various contigs, prep the contigs often from salivary glands, stick them under a microscope and map the location of the probes. Hence, the reluctance to assign a 'genome' based on contigs because the chromosomes cannot be assembled. This is particularly true for a very close relative to the premier model of the Insecta D. melanogaster. A lot of work will have gone into the correct chromosome annotation of this genome. Thus, you can't simply 'promote' a genome without further genomics analysis and it will be a community decision based on the analysis (i.e. reviewers).
Bioinformatic assemblies often get it right. What I would do is use a reference genome assembly based on the nearest ancestor, which you'll know more than me, I'm pretty sure its not D. melanogaster. This is within the context of a comparative genomic analysis between the species/subspecies associated with D. subobscura and then say to NCBI 'thats a genome'. You then have to publish the analysis and submit your assembled genome to NCBI. Its not trivial, but if you 'live next door' to one of the genetic models of eukaryotes (well metazoa) the standards of annotation and assembly will be high.

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


text:
How is an X chromosome encoded into a fasta string?
A human "X" chromosome has a centromere and two "identical" chromatids. If the chromatids are not identical, this fact is not assembled, correct?
The fasta string for a chromomsome is actually representing just one of the chromatids, correct?
So from my understanding, for diploids, the fasta chromosome is actually a mix of 4 chromatids: 2chromatids * 2 pairs of chromosomes into one string. Which for some heterozygous nucleotides (same percentage of each allele) will just be randomly selected from the reads as the correct one.
(unless we are talking about some methods like 10x that can distinguish between chromosomes pairs)
I was thinking which side of the X is read first into the fasta string, and I came to the above conclusion.

So, the two chromatids you see only during meta-phase of replication. In more "regular" state of a cell, there are only the two haplotypes (in the case of a diploid organism)

Well, and regardless of the ploidy, we usually aim for a haploid (monoploid) representation, so it is usually exactly as you say - an arbitrary one of the two states on each heterozygous site is chosen.
However, there are also approaches to phase genomes (to sort which alleles are on which haplotypes). That is possible not just with 10x, but also with long reads (to certain extend) or using trio-binning - a strategy to sequence two parents and the offspring.
However, even phased genomes still contain the reference (fasta) in a haploid state and usually provide the alternative haplotype as a second fasta file.
For the record, these apply to all chromosomes, not just the X.


text:
What visualization was used for ICA in the tissue experiment?
I want to run the experiments in mentioned in enter the matrix.
For the the plot in figure 3B it mentioned that

The relative absolute value of the magnitude of each element in the
pattern matrix indicates the extent to which that sample contributes
to the corresponding source of variation. The sign of the values
indicate over- or underexpression in that factor depending on the sign
of the corresponding gene weights in the amplitude matrix. As a
result, the values can be plotted on the y axis against known
covariates on the x axis to directly interpret the relationship
between samples. When applied to GTEx, we observe one pattern
associated with cerebellum, another pattern that has large positive
values for one donor and large negative values for another donor, and
eight other patterns associated with other sources of variation
(supplemental information online).

I calculated the ICA with head()
        PC1       PC2
0  1.317276 -0.682285
1  1.696110 -0.549004
2  0.037709 -0.548736
3  0.702004 -0.545255
4  1.227363 -0.574219

but fail to see how the authors plotted their ICA graphs


text:
How to add N bases to PyRange object (i.e., shifting PyRanges coordinates)
I would like to add n (e.g., 1k) bases to a PyRange object (i.e., shifting PyRanges coordinates).
How can I do that?
Input:
import pyranges as pr
gr = pr.random(3)

# gr:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10000000 |  20000000 | -            |
# | chr2         |  20000000 |  30000000 | +            |
# | chr22        |  30000000 |  40000000 | -            |
# +--------------+-----------+-----------+--------------+

Desired output:
# magical command[gr + 1000]:
# +--------------+-----------+-----------+--------------+
# | Chromosome   |     Start |       End | Strand       |
# | (category)   |   (int32) |   (int32) | (category)   |
# |--------------+-----------+-----------+--------------|
# | chr1         |  10001000 |  20001000 | -            |
# | chr2         |  20001000 |  30001000 | +            |
# | chr22        |  30001000 |  40001000 | -            |
# +--------------+-----------+-----------+--------------+

You can use the .extend() method. Starting with your data above, gr.extend({"3": -10000, "5": 10000}) would generate the following:
+--------------+-----------+-----------+--------------+
| Chromosome   |     Start |       End | Strand       |
| (category)   |   (int32) |   (int32) | (category)   |
|--------------+-----------+-----------+--------------|
| chr1         |  10010000 |  20010000 | -            |
| chr2         |  19990000 |  29990000 | +            |
| chr22        |  30010000 |  40010000 | -            |
+--------------+-----------+-----------+--------------+

I don't know a way of to apply this "range extension" in a strand specific way (to be able to generate the desired output exactly) but if this is really needed, you should be able to generate 2 gr objects, one for each strand, do the manipulation and bring these together.

OP is not doing anything strand-specific, just adding 1000 to coordinates irrespective of strand. Your answer provides more flexibility than OP needs.


text:
extracting raw gene counts based on the cluster generated by seurat
I'm new to single cell sequencing analysis. I have made different clusters using seurat. Now I would like to extract the raw gene counts based on the cluster generated by seurat.
i have provided my R code used in the analysis.
library(dplyr)
library(Seurat)
library(patchwork)
ftp="https://ftp.ncbi.nlm.nih.gov/geo/series/GSE201nnn/GSE201254/suppl/GSE201254_count_mat_Control-NoCrush.csv.gz"
counts=data.table::fread(ftp)
C_nocrush=counts %>%
select(contains(c("V1","CtCReRGCS")))
rownames(C_nocrush)=C_nocrush$V1
C_nocrush$V1=NULL
CNC=CreateSeuratObject(counts = C_nocrush,project = "He RGC Data", min.cells = 3, min.features = 200)
CNC[["percent.mt"]] = PercentageFeatureSet(CNC, pattern = "^MT-")
VlnPlot(CNC, features = c("nFeature_RNA", "nCount_RNA","percent.mt"), ncol = 3)
plot1 <- FeatureScatter(CNC, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(CNC, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1+plot2
CNC = subset(CNC, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
CNC = NormalizeData(CNC, normalization.method = "LogNormalize", scale.factor = 10000)
CNC = FindVariableFeatures(CNC, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(CNC), 10)
plot1 <- VariableFeaturePlot(CNC)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1+plot2
all.genes <- rownames(CNC)
CNC <- ScaleData(CNC, features = all.genes)
CNC <- RunPCA(CNC, features = VariableFeatures(object = CNC))
VizDimLoadings(CNC, dims = 1:2, reduction = "pca")
DimPlot(CNC, reduction = "pca")
DimHeatmap(CNC, dims = 1:15, cells = 500, balanced = TRUE)
CNC <- JackStraw(CNC, num.replicate = 100)
CNC <- ScoreJackStraw(CNC, dims = 1:20)
JackStrawPlot(CNC, dims = 1:15)
ElbowPlot(CNC)
CNC <- FindNeighbors(CNC, dims = 1:10)
CNC <- FindClusters(CNC, resolution = 0.5)
head(Idents(CNC), 5)
CNC <- RunUMAP(CNC, dims = 1:10)
DimPlot(CNC, reduction = "umap")
saveRDS(CNC, file = "Desktop/Single_cell_RNA/CNC.rds")
cluster2.markers <- FindMarkers(CNC, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)
new.cluster.ids=c("T5/F RGCs", "T5/F RGCs","T-RGCs","Other RGcs and cartpt RGCs","T5/F RGCs","IP-RGCs","Other RGcs and cartpt RGCs","N-RGCs","Alpha-RGCs" )
names(new.cluster.ids)=levels(CNC)
CNC=RenameIdents(CNC,new.cluster.ids)
 DimPlot(CNC,reduction = "umap", label = T, pt.size = 0.5)+NoLegend()

Tried using this command to extract the raw counts based on the cluster but is not working for me.
meta.data.cluster=unique(x=CNC@meta.data$active.ident)
meta.data.cluster=unique(Idents(CNC))
for(group in meta.data.cluster) {
    group.cells <- WhichCells(object = CNC, idents = group)
    data_to_write_out <- as.data.frame(GetAssayData(CNC,slot = 'counts')[,group.cells])
    write.csv(data_to_write_out, row.names = TRUE, file = "1_cluster_outfile.csv")
}

Hi, thank you for supplying detailed code. Could you kindly explain the precise question? I assume the code is throwing an error or not working how it was anticipated. Describing this would be great.

I have updated my script. I want to extract the raw counts of the genes using the selected cluster.

How do you know that the count extraction code you wrote is not working?


text:
Remote blastn - what is breaking my (bash) loop?
I am trying to search the NCBI non-redundant database for sequences similar to a few other (~40) sequences that I already have.
So I've tried running blastn remotely, and looping through multiple queries. I've tried it both as a job (via slurm) and from the command line.
In addition, for every query, I am trying to run 3 searches with slightly different parameters (in WebBlast, I noticed that I recover more sequences of interest this way).
I am trying to dump the output data in a single file (stuff like duplicates can be sorted out later).
This is my code:
for i in /directory/*.fa;do #this is the directory containing the queries
j=${i#"/directory/"}
k=${j%".fa"} #getting a short version of the name for the log file
echo "$k" >> "log.log"
blastn -query "$i" -db "nt" -remote -entrez_query 700:1000[slen] -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
blastn -query "$i" -db "nt" -remote -entrez_query 1000:1400[slen] -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
blastn -query "$i" -db "nt" -remote -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"
done

However, it...doesn't seem to work? Either only the first search for the first query is completed, and then the code does nothing else, or nothing happens. (in a rare case, the first two searches worked) There seems to be no clear distinction between 'it works a tiny bit' and 'it doesn't work at all', which leads me to think it might be something outside the code, e.g. connection issues.
Has anyone faced this before, and would you know how to fix it?
Update
I tried to run a single search as well. As well as adding the waiting period. Right now it doesn't work (same issue, no error message, just sort of appearing to run without results). In the link that Scot provided below, from the NCBI, it was mentioned 'We will move searches of users who submit more than 100 searches in a 24 hour period to a slower queue, or, in extreme cases, will block the requests.' I haven't run >100 searches yesterday (was planning to before I knew this, but that's a different issue).
I think I'll limit the searches as recommended and try to run it from a different computer later on.

Hi @Laura could you have one final try using sleep 60? The other idea is to set up a cron job and run it over night, I wouldn't try and manually overcome the challenge.

@M__ I tried. First search works as normal (finishes in a couple of minutes). The next are...again, very very slow. I left the job running for a bit, and it's taken 2h for now while only on the 4th search.

Hi @Laura - depending on your network environment (perhaps on a University network), and the details of how NCBI BLAST throttling works, you may be throttled as a group rather than as an individual. Trying this from a different network may help if this is the case.

And it is beyond the scope of your question, but if you are desiring to run hundreds of BLASTs eventually, I'd look into downloading the BLAST executable and the nt database, and running the alignments locally. nt is not a small database, so this may be best done using slurm on your cluster, but if you can limit the alignments to a taxonomic segment of nt, you may be able to do this locally using just that taxonomic slice of nt.

Answered in the post

I gave it a whirl using a random bit of SARS-CoV-2
blastn -query out.fa -db "nt" -remote  -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results.tsv"

blastn -query out.fa -db "nt" -remote -perc_identity 94 -outfmt "6 sacc sseqid sseq evalue bitscore length" >> "results2.tsv"

Both gave a 500 line output - the limit of the blast specified.
One possibility is to place a sleep 10 line between each blastn command (3 additional lines).
The rationale is specifications NCBI set:

Do not contact the server more often than once every 10 seconds. Do
not poll for any single RID more often than once a minute. Use the URL
parameter email and tool, so that the NCBI can contact you if there is
a problem. Run scripts weekends or between 9 pm and 5 am Eastern time
on weekdays if more than 50 searches will be submitted.

The suggestions haven't worked, so I'd recommend trying sleep 60 between each blastncommand on the basis this might be an issue with RID. If that fails the fall back is to run the job overnight US time via cron tab, which I'll post the code if it is required. Thus,

Run scripts weekends or between 9 pm and 5 am Eastern time
on weekdays if more than 50 searches will be submitted.

The final thing to add is that it is a fairly long job so connectivity might be an issue. If this was the case then it might also be resolved by running the script overnight in the OPs timezone, i.e. via cron

sleep 60 works but is very slow.

Agreed one option is to download the database and perform the Blast locally.
Alternative use mmseqs2
You could perform the Blast on the Cloud, like GCP (much faster).
Personally I would look at performing over night runs when the network was not in use.

FYI @scot posted first ... similar ideas but credit to them.

This is likely a throttling issue - see here for some details, specifically the one about contacting the NCBI server more than once every 10 seconds.
You may be able to submit a multiFASTA of all your queries to get around this issue, rather than individual FASTA sequences, or alternatively, add in a sleep 10 command between your BLAST commands and at the end of each for loop iteration to throttle the requests.


text:
Efficiently count reads overlapping a split feature
I have a bam file with reads, and a set of features. The particularity is that some features are "split", being made of several sub-features. For example, here bed1 and bed2 are made of a simple feature, bed3 contains a split feature:
     0        10        20
     |---------|---------|
bam   >>
       >>>>>>>>>>>>
           >>
     |---------|---------|
bed1 ####
bed2             #####
bed3 ####        #####

I want to count the reads overlapping with each (set of) feature, but importantly I want to count only once a read that intersects both features in a set:
feature   count
   bed1       2
   bed2       1
   bed3       2   -> NOT 3

In R, I can do that with a GRangesList where bed1 and bed2 are list elements containing 1 GRange, whereas bed3 is a list element containing 2 GRanges; find the reads that intersect a list element, then use unique() to avoid double-counting reads that intersect two GRanges in the same list element:
suppressPackageStartupMessages(library(GenomicAlignments))
exple_beds <- GRangesList(list(GRanges("I",
                         ranges =  IRanges(start = 0, end = 3),
                         strand = "+",
                         names = "bed1-a"),
                 GRanges("I",
                         ranges = IRanges(start = 12, end = 16),
                         strand = "+",
                         names = "bed2-a"),
                 GRanges("I",
                         ranges = IRanges(start = c(1,12), end = c(3,16)),
                         strand = "+",
                         names = c("bed3-a", "bed3-b")))) |>
  setNames(c("bed1","bed2","bed3"))

exple_bam <- GAlignments(seqnames = Rle(rep("I", 3)),
            strand = Rle(rep(factor("+", levels = c("+","-","*")),3)),
            pos = c(1L, 2L, 6L),
            cigar = c("2M","12M","2M"))

sapply(exple_beds,
       \(.x) findOverlaps(.x, exple_bam) |>
         to() |>
         unique() |>
         length())
#> bed1 bed2 bed3
#>    2    1    2

Created on 2022-10-07 by the reprex package (v2.0.1)
However this gets really slow with a big BAM file.
So, is there

either a more efficient way to do it in R,
or a way to do it with bedtools and friends without creating 100,000 BED files?

For anyone reading this later, I did not find a satisfying answer and ran the above R code on subsampled bams to allow it to run in a reasonable duration (in my case it's "good enough"). An easy improvement is to parallelize the sapply(). For Python aficionados, I expect RSeQC is doing something similar.


text:
What can be the bias of aligning paired-reads in a single-end mode?
I don't know if I'm in the right place but I have a technical problem to fix. I would have to align paired-end reads from Illumina sequencing to compare a normal genome with a tumor one.
When I align with bwa-mem (default parameters) with paired-end mode I get a bigWig file where the coverage of the two samples is quite similar. When I align with bwa-mem (default parameters) with single-end mode I get coverage holes (definitely does not align any read) in some specific regions. But this only happens in tumor genomes.
Is it possible that something is really going on in that region? Or is it just a bias? And if it's a bias, why in the normal condition, in that specific region, do they align?
I think I understand that by treating paired-end reads as single-end I lose the directionality of mate2 if I align mate 1. But I still can't understand why I have this alignment defect only in the tumor genome and not in the normal genome (I observe this phenomenon in all the samples analyzed, not just in one).
Could someone give me some explanation or recommend any checks to do?
spoiler - these holes fall exactly in highly repetitive regions

In principle one should simply not align paired-end reads as if they were single end reads, not only because you are not taking advantage of the forward/reverse information contained in the reads, but also you are not taking into account the expected distance between the two pairs, both of which result in better mapping rates/accuracy.
Having said that, and not recalling from the top of my head right now the relevant BWA-mem algorithmic implementation details, it is generally difficult to map (even) paired-reads to repetitive regions. Still there will be some rate of success. In tumor samples however, we could reasonably imagine that chromosomal rearrangements, INDELs, CNVs, etc. will make it even harder to confidently assign these reads to a location in the reference, since these extra variation will increase e.g. multiple mapping locations and mapping qualities for some fraction of the reads, thus resulting in regions with low or no coverage compared to non-tumor samples.
Repetitive regions might exacerbate the problem and get affected even further compared to other regions.

I tried to give me this explanation: with single-end alignment of paired-end reads, losing the directionality of the reads and also the expected fragment length, the R2 no longer maps if it maps the R1, but can map as a single reads anywhere, or don't map at all.
I thought that therefore I don't see alignment (in tumors) in repetitive regions because of the high mutation rate I have in the tumor sample. While in normal I don't have this situation.
However, aligning with paired-end reads, I see coverage in these regions again despite the high mutation rate because R1 and R2 help each other map in the best position, and so even if I have mutations I still see the alignment.
Could it be appropriate to make a variant call despite being on repetitive regions?

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.

Check on the topic of 'mapability' (e.g. 1000 genomes or UCSC tracks, snpable, ..). Often it is practice to mask out regions of poor mapability in order to decrease the rate of false positives. Repetitive regions fall often within poorly mappable regions, your read length becomes somewhat relevant here


text:
Name this alternative INDEL format
I have recently been working with variants which are represented in a non-standard format. This alternative format uses a period instead of prefixing variants with a reference base. What is this format called?
There is no specific output format associated with it. It's part of some calling code I've inherited. I figured that this existed elsewhere. But, I could be wrong.
Expected Format:
Insertion:
pos=4
alt=ATTT
ref=A

Deletion:
pos=4
alt=A
ref=ATTT

Alternative Format:
Insertion:
pos=5
alt=TTT
ref=.

Deletion:
pos=5
alt=.
ref=TTT

The alternative deletion representation that you are describing, i.e. using . to represent alternate alleles that are not present, is allowed within the definition of the VCF format (in other words, this representation can be found in a standard format). An example of this is provided on the Wikipedia page for VCF:
#CHROM POS      ID         REF   ALT    QUAL  FILTER
20     1230237  .          T     .      47    PASS

While missing values are allowed for ALT, they are not typically used because they lead to ambiguous variant definitions.
For missing REF values, the VCF format specification has the following explanation:

Each base must be one of A,C,G,T,N (case insensitive). Multiple bases are permitted. The value in the POS field refers to the position of the first base in the String. For simple insertions and deletions in which either the REF or one of the ALT alleles would otherwise be null/empty, the REF and ALT Strings must include the base before the event (which must be reflected in the POS field), unless the event
occurs at position 1 on the contig in which case it must include the base after the event....

To work around this ambiguity, VCF files are commonly left-normalised so that the specified sequences and variants can be precisely placed within the reference genome. Left normalisation can be carried out using bcftools norm.
More details on left normalisation can be found here:
https://genome.sph.umich.edu/wiki/Variant_Normalization#Left_alignment

A variant is normalized if and only if it is parsimonious and left aligned.
A variant is parsimonious if and only if it is represented in as few nucleotides as possible
without an allele of length 0.
A variant is left aligned if and only if it is no longer possible to shift its position to the left while keeping the length of all its alleles constant.

In other words, these are non-parsimonious variants.


text:
Why does the Smith-Waterman algorithm have a complexity of $O(m^2 \times n)$?
I recently read that original the Smith-Waterman algorithm to do a local pairwise alignment of sequences has a complexity of $O(m^2 \times n)$, where $m$ and $n$ would be the lengths of the sequences with $m$ the potentially longer one. What I am wondering is where the additional complexity arises in the algorithm compared to the Needelmann-Wunsch algorithm of global alignment, which has complexity $O(m \times n)$.
Am I correct to assume the $m\times n$ for the Needelmann-Wunsch algorithm comes from building up the scoring matrix - whereas the traceback is negligible in complexity i.e. $O(m\times n + m)$? Hence we arrive at $O(m\times n)$.
Why is there a $m^2$ term for the Smith-Watermann algorithm?

I think your source is wrong. You are correct to assume the brunt of its complexity comes from building the $m *n$ matrix. The main differences are in the calculations per cell and the traceback, but these don't make $m ^2$. Wikipedia also gives an identical $O(m * n)$ as worst case for both algorithms.
edit:
Aparently in this paper the algorithms were improved to go from $m^2n$ to $m*n$

Yes, the Smith-Waterman algorithm is often called the Smith-Waterman-Gotoh algorithm for this reason.

Thank you for the great reference. I gave it a read. So under certain gap penalties, the complexity is reduced. I don't yet fully understand why though. Is there a heuristic explanation?


text:
How to plot stacked bar chart using R showing mean with range and labelled values?
I am trying to plot a stacked bar graph with mean values from multiple repeats, also showing the range of the data and values labelled.
My data frame is below:
structure(list(Average.DB.per.FA = c(0, 0.5, 1, 1.5, 2, 0, 0.5,
1, 1.5, 2, 0, 0.5, 1, 1.5, 2, 0, 0.5, 1, 1.5, 2, 0, 0.5, 1, 1.5,
2, 0, 0.5, 1, 1.5, 2), condition = c("CDD", "CDD", "CDD", "CDD",
"CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD", "CDD",
"CDD", "CDD", "STD", "STD", "STD", "STD", "STD", "STD", "STD",
"STD", "STD", "STD", "STD", "STD", "STD", "STD", "STD"), percentage = c(0,
0.807218011, 0.192781989, 0, 0, 0, 0.737968015, 0.262031985,
0, 0, 0, 0.739096101, 0.260903899, 0, 0, 0, 0.466751289, 0.346156739,
0.170617491, 0.016474482, 0, 0.462997344, 0.331833715, 0.183359127,
0.021809813, 0.012053568, 0.47623097, 0.307175237, 0.182459607,
0.022080618)), row.names = c(NA, -30L), class = "data.frame")

I tried to use the following code:
    ggplot(long_DF, aes(fill=Average.DB.per.FA, y=percentage, x=condition)) +
  geom_bar(position=position_dodge(), stat="identity")

But the graph I obtained does not have bars add up to 1:

I would like to have graphs like below (which I created by Prism using the same set of data but Prism does not allow me to label bars with numbers):

Is there any way to reproduce the graph by Prism using R and add values to each stacked bar?
Appreciate your help and thanks in advance.

Hi its a great question. Could you kindly explain the biological aspect of the code, i.e. why is this bioformatics rather than pure statistics?

It seems like this is missing a grouping variable in the aesthetics definition. The R help page for position_dodge mentions the following:

Position_dodge() requires the grouping variable to be be specified in the global or geom_* layer. Unlike position_dodge(), position_dodge2() works without a grouping variable in a layer. position_dodge2() works with bars and rectangles, but is particulary useful for arranging box plots, which can have variable widths.

Bearing this in mind, setting group=Average.DB.per.FA seems to give the result you want:
long_DF %>%
  ggplot() +
  aes(fill=Average.DB.per.FA,
      y=percentage, x=condition, group=Average.DB.per.FA) +
  geom_col()

... although this is adding up to 3 rather than 1, as would be expected from looking at the category counts:
> long_DF %>%
+   group_by(condition) %>%
+   summarise(sumPct = sum(percentage))
# A tibble: 2 × 2
  condition sumPct
  <chr>      <dbl>
1 CDD            3
2 STD            3

Note: as seen above, I'm using geom_colinstead of geom_bar with a different statistical summary method. If you're using stat=%something%, there's probably a better way to approach the task. The following is from the first paragraph of the R help page for geom_bar/geom_col:

There are two types of bar charts: geom_bar() and geom_col(). geom_bar() makes the height of the bar proportional to the number of cases in each group (or if the weight aesthetic is supplied, the sum of the weights). If you want the heights of the bars to represent values in the data, use geom_col() instead. geom_bar() uses stat_count() by default: it counts the number of cases at each x position. geom_col() uses stat_identity(): it leaves the data as is.


text:
Running a Python Script from a bash script on multiple input files across different directory paths
I am in the process of writing up a bash script that runs three bioinformatics tools and then executes a python script. This python script compiles some of the data from the informatics tools output into a single merged output file. Now I need to run this on a large number of samples at a time so I am trying to set up the script to read through the directory path to where all the output from the tools is stored.
These are the directory paths for the python script input files:
EFS/Analyses/BGC.test/Sample_Name/antismash2/sample.LargeContigs.gbk
EFS/Analyses/BGC.test/Sample_Name/Bagel4/00.OverviewGeneTables.json
EFS/Analyses/BGC.test/Sample_Name/Deepbgc/Deepbgc.bgc.tsv

Now there are multiple Sample_Name directories containing identical structures in the EFS/Analyses/BGC.test path. I would like to run the script on all of them
To do this, I have written a for loop within the bash script as follows:
sample_paths=$(find $EFS/Analyses/BGC.test -mindepth 1 -maxdepth 1 -type d \ )
while read sample_paths;
do
   sample=$(basename $sample_paths)
   GenBank_Input="$sample/antismash2/$sample.LargeContigs.gbk";
   TsvFile_Input="$sample/Deepbgc/Deepbgc.bgc.tsv";
   Json_Input="$sample/Bagel4/00.OverviewGeneTables.json";
done <<< “$sample_paths”

export GenBank_Input
export TsvFile_Input
export Json_Input

The goal here being to define the sample path for every sample and then for every sample, looping through to the specific directories of each input file. I then exported the variables and used os.environ from the os package to use them in the python script as such:
    #parse gbk input file
    for record in SeqIO.parse(os.environ["GenBank_Input"], 'genbank'):
    #parse tsv input file
    file_handle = open(os.environ["TsvFile_Input"])
    #parse json input file
    json_data = json.load(open(os.environ["Json_Input"]))

Now the error I get when I run the script is a FileNotFoundError: [Errno 2] No such file or directory: '“”/antismash2/“”.LargeContigs.gbk' so the error seems pretty clear cut. The issue seems to be with the path that I defined within the while loop variables but as of right now, I have been unable to resolve. I am trying to familiarize myself more with bash scripting as I am still fairly new to it. Any guidance on how I might adjust the code to resolve this and get the script working is much appreciated, thanks!

Could you kindly print the output of sample_paths @Joe_G?

What is $EFS? How and where is that defined? Also, note that you have “$sample_paths” (curly, UTF quotes) instead of "$sample_paths" (normal ASCII quotes). Finally, are you expecting these variables to be strings or arrays? Bash strings interpreted as python lists? Something else?

Your loop is more complicated than you need and it won't work properly anyway since you're using the same variable (sample_paths) to iterate and save. Try this instead:
for dir in EFS/Analyses/BGC.test/*/; do ... done

Next, all you really need is the sample names. The rest can be inferred:
sample_names=""
for dir in EFS/Analyses/BGC.test/*/; do
    sample_names+="$(basename "$dir") "
done
## remove the extra space added in the end
sample_names=$(sed 's/  *$//' <<< "$sample_names")
export sample_paths

Then, in python, you would do:
#!/usr/bin/env python3

import os

base_dir="EFS/Analyses/BGC.test"
for sample_name in os.environ["sample_names"].split(' '):
    gbk_file = """%s/%s/antismash2/%s.LargeContigs.gbk""" % (
        base_dir, sample_name, sample_name
    )
    tsv_input = """%s/%s/Deepbgc/Deepbgc.bgc.tsv""" % (base_dir, sample_name)
    json_input = """%s/%s//Bagel4/00.OverviewGeneTables.json""" % (base_dir, sample_name)
    print("============== Sample Name:",sample_name, "=============")
    print("GenBank file:", gbk_file)
    print("Tsv file:", tsv_input)
    print("Tsv file:", json_input)

I created the following test directories:
$ tree
.
└── EFS
    └── Analyses
        └── BGC.test
            ├── Sample_1
            │   ├── antismash2
            │   │   └── sample.LargeContigs.gbk
            │   ├── Bagel4
            │   │   └── 00.OverviewGeneTables.json
            │   └── Deepbgc
            │       └── Deepbgc.bgc.tsv
            ├── Sample_2
            │   ├── antismash2
            │   │   └── sample.LargeContigs.gbk
            │   ├── Bagel4
            │   │   └── 00.OverviewGeneTables.json
            │   └── Deepbgc
            │       └── Deepbgc.bgc.tsv
            └── Sample_3
                ├── antismash2
                │   └── sample.LargeContigs.gbk
                ├── Bagel4
                │   └── 00.OverviewGeneTables.json
                └── Deepbgc
                    └── Deepbgc.bgc.tsv

And running the two snippets of code above produces:
$ foo.py
============== Sample Name: Sample_1 =============
GenBank file: EFS/Analyses/BGC.test/Sample_1/antismash2/Sample_1.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_1/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_1//Bagel4/00.OverviewGeneTables.json
============== Sample Name: Sample_2 =============
GenBank file: EFS/Analyses/BGC.test/Sample_2/antismash2/Sample_2.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_2/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_2//Bagel4/00.OverviewGeneTables.json
============== Sample Name: Sample_3 =============
GenBank file: EFS/Analyses/BGC.test/Sample_3/antismash2/Sample_3.LargeContigs.gbk
Tsv file: EFS/Analyses/BGC.test/Sample_3/Deepbgc/Deepbgc.bgc.tsv
Tsv file: EFS/Analyses/BGC.test/Sample_3//Bagel4/00.OverviewGeneTables.json

All this said, I don't really understand why you are using the shell at all. You can implement all of this directly in python:
#!/usr/bin/env python3

import os

base_dir="EFS/Analyses/BGC.test"
sample_names = os.listdir('EFS/Analyses/BGC.test/')

for sample_name in sample_names:
    gbk_file = """%s/%s/antismash2/%s.LargeContigs.gbk""" % (
        base_dir, sample_name, sample_name
    )
    tsv_input = """%s/%s/Deepbgc/Deepbgc.bgc.tsv""" % (base_dir, sample_name)
    json_input = """%s/%s//Bagel4/00.OverviewGeneTables.json""" % (base_dir, sample_name)

    # do something with these files here

thanks for the reply. You're right, a direct python implementation would be much easier. I was trying to familiarize myself some more with utilizing bash as well. That being said, I will try implementing the python loop you provided. Now if I have extraneous files in the main BGC.test directory, I will have to ignore these manually?

@Joe_G no, only if you have extraneous directories. The trailing slash in for dir in EFS/Analyses/BGC.test/*/ limits the results to directories only. You can do the whole thing in bash too, it's the mixing of the two languages that makes it more complicated.

Gotcha, so If I were implementing directly through python (ie. negating bash all together), does that sample_names = os.listdir ... you provided in that example python code accomplish essentially the same thing? (ignoring files and looking just at directory paths)

@Joe_G no, that doesn't ignore files. But you can just test the results and skip any that are files or put the whole thing inside an if os.isdir(sample_name_):.


text:
Regarding chi square?
I need some suggestions regarding an analysis. I might be wrong as I have no stats background.
I have list of exon with differential usage across samples. I can find lots of denovo mutations from both proband and siblings in those exons. Now I am trying to check was there any enrichment of proband specific denovo mutations in those exons. I can think of doing a chi square test. But I am not sure what will be expected value in chi square test. Imagine I am having 700 denovo proband mutation in 650 samples out of 1500 probands and 300 denovo mutations in 320 unaffected  siblings, out of 1200. Any suggestions? Thanks

Denovo
WT
Total

Proband
553
1202
1755

Unaffected sibs
505
813
1318

Total
1058
2015
3073

Hi @user3377241, I tried to convert your numbers into a table, but it seems inconsistent / incomplete. Could you please fill in the table? I think this will help us to understand what's going on.

@gringer I have edited the table with actual numbers. These number are of probands and siblings. The variants number will be different as many a time a proband have multiple events. I am not sure with should be expected value in these calculation. Thanks

I'm a bioinformatician, but not a statistician, and it's been a while since I did any $\chi^2$ tests, so take this answer with a grain of salt.
From what I remember, the most important thing is to be clear about what your observed and expected values are, and to be consistent. In other words, it doesn't matter what you pick, as long as you make it clear what you've picked. For your situation there are at least four possible "expected" values that can be used:

Average counts based on your available data (this is probably the most commonly used approach)
Equal proportions based on your available data (this is less common, because it alters table totals)
Unaffected sib counts (also less common, because it alters totals)
Population average counts based on public data

Unfortunately, there's no answer that will work for all situations. It's a good idea to think about which of these (or any others) makes sense in your particular case, what you're actually trying to compare, and how much you want to generalise to other populations. It's almost always the case that you know your own data better than anybody else, so are in the best position to make this decision. For example, if this is a very rare variant in the general population, you might choose to look for general population data for the variant, but if you're only interested in the impact of people who are aware of a family history, then unaffected sibs might be a better approach.
In any case, the calculations for generating the expected values are similar: work out the proportion of the variant in the expected population, and apply that proportion to the total count in the observed population.
In the first case, the expected counts would be something like this [note that the table totals are preserved]:

Denovo
WT
Total

Proband
*604.2
1150.8
1755

Unaffected sibs
453.8
864.2
1318

Total
1058
2015
3073

* calculation: row total (1755) as a proportion of full total (3073) = 0.5711, which is then multiplied by the column total (1058) = 604.2
In the second case, it would be this [note that the row totals are different]:

Denovo
WT
Total

Proband
529
1007.5
1536.5

Unaffected sibs
529
1007.5
1536.5

Total
1058
2015
3073

In the third case, it's simply comparing Proband vs Unaffected counts, normalising based on the row total:

Denovo
WT
Total

Proband
672.4
1082.6
1755

Unaffected sibs
505
813
1318

Total
1177.4
1895.6
3073

Note that the second row can be omitted from $\chi^2$ calculations here, because it's identical to the expected values, so adding that difference will not change the calculated $\chi^2$ sum.
[this could also be done based on the column total, if you wanted to compare Denovo vs WT, instead of Proband vs Unaffected]

It seems reasonable, it depends how you do the Chi2.
The Chi2 distribution depends on a good estimate of the expected (E). Here your observed (O) is the "affecteds" and you want to assess whether this has a suspected elevated  mutation rate, then that looks a meaningful measure.
There's two approaches

Your E is the unaffected siblings and E is calculated by averaging this value. You could also perform the reverse calculation where E was the average for the affecteds and O is the unaffected. There will be an issue about the degrees of freedom, which would require some thought.
There is an alternative way to perform Chi2, which I forget because O and E are matched (both siblings), here there is no averaging. You'd need to look it up, the numerator is the summed difference for each sibling group ... but I can't remember.

Point 2 would be a good test, but you'd need to check the test.
A phylogeneticist would also use the Chi2 in a very different test. They would draw a tree and designate two rates, one for affected and unaffected. If there was elevated mutation rate in the affected it could be tested on a likelihood ratio test, i.e. the rate of evolution was increased. This can be done because the test tree would be one where the rate was the same between affecteds and unaffected progeny. The way the test works is the likelihood difference between the two trees is assessed and because the number of increased parameters (degrees of freedom) is known and log probabilities can be subtracted from each other and directly assessed (I think its degrees of freedom + 1 ... 1 additional parameter so its 2) - and therefore the difference in likelihood is tested directly against the Chi2. Its kinda neat and exploits the behaviour of log probabilities.


text:
Is there an equivalent of Entrez queries for local BLAST?
I have been having some trouble with running command line BLASTn remotely. Therefore, upon recommendations, and considering I have >100 searches planned, I am considering just downloading the nt database (thankfully possible).
My issue is that one of the settings that proved valuable with my data in WebBlast is running an Entrez query. I get more useful results for each query when instead of one search, I run three, with each of the sub-searches delimited by subject sequence length (e.g. -entrez_query 700:1000[slen] to only get the results which are between 700 and 1000 nucleotides long).
WebBlast and remote Blast have that, local Blast does not.
Since I only really need 2 Entrez queries, I guess it would not be conceptually impossible to get my copy of nt database as FASTA, and then run a script to split it into three directories (700-1000nt, 1000-1400nt, the rest while discarding all<700 nt), then make three blastable databases out of it. However it does feel rather inelegant. I was wondering if anyone might have a better solution!

Comments are not for extended discussion; this conversation has been moved to chat.

I looked back at the question and have a more comprehensive solution.

There is no option for -entrez_query on a local blast search. The information is here

entrez_query   Restrict search with the given Entrez query. Remote searches only.

There is an option without -outfmt for length and that is worth investigating.

length means Alignment length

Perl one liner solution to NCBI's many parameter problem

a) Preprocess your fasta file into separate databases of the required length ... by turning a fasta flatfile db into one where all the subject sequences are 700 - 1000.
perl -0777 -nle 'while (/^(>.+\n^[AGCT-NRYMKWS]{700,1000})\n/gm){print $1,"\n"}' example.fa > example700-1000.fa

b) Now make the database on the restricted length fasta flatfile db,
makeblastdb -in example700-1000.fa -parse_seqids -blastdb_version 5  -out laura700-1000.db -title "Laura 700 - 1000" -dbtype nucl

This db will now contain an equivalent of the -entrez_query 700:1000[slen]. Obviously, you would need to make a new db for each change in the sequence length query. Note the option -taxid_map.
c) Finally, run the blast search on the restricted length blast db.
blastn -query laura.fas -db laura700-1000.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

Code testing
Input is test.fa
>here
AGAGACGTAGAGAAAAAAAAAA
>here2
TAGAGA
>here3
AGAGACGTAGAGA

Modified perl oneliner because sequence length is different,
perl -0777 -nle 'while (/^(>.+\n^[AGC-T]{6,13})\n/gm){print $1,"\n"}' test.fa

Output
>here2
TAGAGA
>here3
AGAGACGTAGAGA

Thus the sequence that was outside the range 6 to 13 (sequence here2 and here3) is excluded (here). This seems alot easier to me that all those blastn options.
Note NRYMKWS refer to the degenerate nucleotide code

"Often, one needs to search ... a specific subset of sequences within an existing database. For these type of searches a convenient way to conduct them is by creating a virtual BLAST database"
From "Use blastdb_aliastool to manage the BLAST databases" https://www.ncbi.nlm.nih.gov/books/NBK569848/


text:
Summing up data effectively in excel
I have some  tsRNA count data in excel which I need to add selectively in blocks based on the count
refID   Count
1) tRNA-Ala-AGC-1   1
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  1
tRNA-Ala-AGC-1  3
tRNA-Ala-AGC-1  300
tRNA-Ala-AGC-1  51
tRNA-Ala-AGC-1  86
tRNA-Ala-AGC-1  6
tRNA-Ala-AGC-1  2
tRNA-Ala-AGC-1  1
tRNA-Ala-AGC-1  5

2)tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 2
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 64
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 2
tRNA-Ala-AGC-11 3
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1
tRNA-Ala-AGC-11 1

3) tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3
tRNA-Ala-AGC-3

I am using sum and subtotal functions in excel to get the total counts of each of the blocks.
I have about 22,000 values and several hundred groups like this.
Is there a way to sum up all such blocks all at once based on the groups or at least a shorter way to sum things work other than autosum and subtotal?
Any help will be useful.

Please don't do bioinformatics in excel.

Thanks for including the input here, but it doesn't look like what I'd expect to see in a well-formatted Excel spreadsheet. Are those X) really in the file? Why are there no counts for tRNA-Ala-AGC-3? Does your input table have empty lines?

Suppose the header is on row 1 and starts at column 'A'. Then, for the above example, just point to a cell at say C1 to C3 and write ...
=sum(B2:B12)
=sum(B14:B25)
=sum(B27:B38)

This gives the subtotals.
For the total sum point the mouse at D1 and write:
=sum(B2:B38)

... except of course the real value is ...
=sum(B2:B22000)

In all cases you must put the = to start the cell.
Its better to try out R, or Python pandas or just a Python list print (sum([1,2,3,4,5]) but Excel's okay.

The issue is its slow, its taking hours to calculate the sum of 22 000 cells - thats not a lot even for Excel. In Python that wouldn't even exceed 1 second, literally. What you are facing is a RAM bottleneck - there is literally no other explanation.
Two suggestions firstly, look at Windows Excel, particularly if you are using OSX. This is because historically Windows was much more RAM efficient than OSX. Now however, I don't know. Secondly, simply switch to a desktop with 4 Gig RAM or more. Literally thats all it will take.
The drop back position is to do the calculation in two halves on the low RAM machine, i.e. split the data set in two. Not a cool solution but it will probably work. The final solution is to look at the stuff below the line.

The solutions are to shift to VB which none of us code in these days, but its the code underneath Excel, or the alternative is to write a custom script in a different language that will do the adding such as Python, R ... even bash.

I have already done that but since there are so many cells I have to do the sums as you suggested for each subgroup individually and it takes several hours to do such sums, hence I am looking for a quick work around.

Hi @AranyakGoswami what you are up against is a RAM bottleneck. Excel shouldn't have trouble with the calculation on a standard desktop.

I suspect @AranyakGoswami is saying that it takes several hours to select regions around the counts for each subgroup, and add the associated formula into the correct cell near the gene counts. At about 1 minute per gene, several hours could mean about 10,000 genes.

It looks like what you want to do is a pivot table, which provides summary statistics for fields, grouped by specified categories. As a prerequisite, this requires an input dataset that is well formatted; more generally:

Put each observation or sample in its own row; don't leave row gaps between samples or groups. Leaving gaps makes bulk analysis more difficult and time consuming. It can still be done, it just means that the data needs to be pre-cleaned, or alternatively analysed.
Put all your variables in columns - the thing that vary between samples, like ‘strain’ or ‘DNA-concentration’
Have column names be explanatory, but without spaces
Do not combine multiple pieces of information in one cell

For more information, see this Data Carpentries Workshop on the recommended way to structure data in spreadsheets:
https://datacarpentry.org/organization-genomics/01-tidiness/index.html#structuring-data-in-spreadsheets
If the data is prepared in this way, then selecting the whole table and creating a pivot table should work. Unfortunately, the programs keep changing where in the menu that can be found. Currently in LibreOffice, it's Data -> Pivot Table -> Insert or Edit..., wheareas on my version of Excel on Mac OSX, it's found at Data -> Summarise With Pivot Table.
This will load up a visual interface that allows you to drag and drop fields to summarise. For this data, you'd probably want RefID in rows, and sum(Count) as the data field:

After setting this up and clicking Ok, you should be greeted with a smaller table containing summarised variables:

I have purposely added a blank space after each block so that I may avoid errors for all the sums made.

It may be that the pivot table will work with gaps, but bear in mind that leaving gaps in general makes bulk analysis more difficult. It can still be done, it just means that the data needs to be pre-cleaned, or alternatively analysed. I've added this clarification to my answer.


text:
How to get the external (cluster) job id and not internal snakemake job id in cluster log files?
Background I am using snakemake and I want to add the external (cluster) jobid to my cluster output and error log files. Below is my shell script.
snakemake --configfiles config/config_resources.yml config/config_parameters.yml --latency-wait 60 --use-conda --printshellcmds --jobname PipeName_{name}.jobid{jobid} --jobs 300 --cluster "bsub -W 480 -n {threads} -M {resources.mem_gb}G -R 'rusage[mem={resources.mem_gb}G]' -o logs/{name}_{jobid}.out -e logs/{name}_{jobid}.err".
But I only get the internal one in de log name from snakemake (e.g. 1 to 5), instead of the one from the cluster.
Question How to I get the external job id in there and not the local one made by snakemake? Or where in my shell script did I go wrong?
Cheers and thanks in advance.

bsub supports some filename placeholders using a percent sign, like %J for the LSF job ID, so you can use that:
snakemake ... --cluster "bsub ... -o logs/{name}_%J.out -e logs/{name}_%J.err"
That way snakemake will format the cluster command string with its wildcards (like name) as it constructs the bsub command, and then bsub will put the job ID in place in those filenames.  Since Snakemake just sees that as one big command to run for the cluster and doesn't know the details other than formatting the text, I think the two shouldn't clash.

Oh my, thank you! That worked!


text:
Is it possible to integrate multiple gene expression datasets and use it for WGCNA?
I have 8 RNA-seq datasets and am interested in looking at genes co-expressed with a specific gene.
Among 8 RNA-seq datasets, 6 have less than 20 samples. Rather than working on each dataset individually for WGCNA, do you think I can merge all datasets and use them for WGCNA?

Yes, if the expression datasets are from the same platform, you can merge all the datasets for common genes among them. However, if they are from different platforms then some batch correction needs to be done before analysis. You can plot a PCA to check how much variability there is among datasets from different platforms and then decide whether batch correction is needed.

thanks a lot nitesh

I would try to process all data sets in a similar way, then combine them and batch-correct, e.g., using ComBat or WGCNA's own empiricalBayesLM.

thanks a lot, peter. Can I use edgeR removebatcheffect?

It should work as well.

Hi @Peter could you please help me with this https://bioinformatics.stackexchange.com/questions/19922/how-to-find-closely-related-genes-for-a-specific-gene-from-wgcna-modules


text:
Compare VCF files
I have two strains, X and Y, and I would like to compare the variants found in both. I have my own parent strain, Z, that was sequenced and assembled to be used as the reference genome.
I first aligned X and Y to Z using minimap2 and then I used FreeBayes as the variant caller. Once I had the two separate VCF files I looked online to see what was the best way to filter out the unique and common variants in X and Y.
I found the bcftools isec tool to be my best bet so I went ahead and ran this line of code:
./bcftools isec -p out/path/ X.vcf.gz Y.vcf.gz
Now I have three new VCF files:

0000.vcf (variants unique to X)
0001.vcf (variants unique to Y)
0002.vcf (variants in both X and Y)

However, it's still in VCF format. I need a way to parse through the VCF file to list out the SNPs with high quality with their position on the genome. I tried using bcftools stats but it just lists out how many SNPs there are and not their position.
TLDR: Can someone help me find a way to list out the locations of SNPs from a VCF file in human readable format?

Can you explain what you need a bit more clearly? Try showing us a few lines of the VCF file and what you would like to convert them to, for example. Since the VCF already gives locations of variants in a human readable format, your question is a bit unclear. Do you want SNPs only, without any larger variants? Do you just want the first, second, fourth and fifth fields of the VCF? Do you want to convert it to hgvs notation? Something else?

If I understood you correctly you are searching for a way to represent the genomic locations in your vcf file. However, the vcf file is already a human readable file (vcf specification).
The first lines are the header, but afterwards you can find your information in a tab delimited fashion:
#CHROM POS ID REF ALT QUAL FILTER INFO

So the information on genomic positions is directly readable from the vcf file (first 2 columns)


text:
Stringtie/DEprep.py gene/transcript IDs are wrongly formatted
Hi my RNAseq workflow is ending up with wrongly formatted gene IDs (and separately transcript _IDs) after a Hisat2 ->samtools sort -> stringtie -> DEprep.py workflow.
DEprep.py outputs a gene_counts.csv which includes this...
gene:ENSRNOG00000000017         47
ENSRNOG00000000021              220

Spelling this out "gene:" should not be in the gene_ID.  First thought of course was to just delete the "gene:" string with search and replace.
However, sometimes, the same gene appears as both indicating the error could be really damaging:
Gene_ID                   Counts
ENSRNOG00000004937        5234
gene:ENSRNOG00000004937   0

Digging around I have found that the issue is present in the initial stringtie /ballgown/ GTF output prior to running DEprep.py
chr1    ensembl exon    la la   gene_id
 "gene:ENSRNOG00000046319" la la la

Curiously my workflow is copied from my own local Galaxy installation workflow that uses the same steps together with the SAME fastq's the same, Hisat2 index, the same Guide.GFF, but there is no issue there.
I have spent a day googling this and no success yet, lots of peep have missing ids. One other person had this symptom, but sadly no resolution.  Any suggestions gratefully received!!
Code details below and an attached folder in the comments with dummy data and r script that reproduces this error.
HiSat2:

hisat2 --dta --add-chrname --time --summary-file YH2minihisat.txt
--rna-strandness FR  --pen-cansplice 0 --pen-noncansplice 12 --pen-canintronlen G,-8.0,1.0 --pen-noncanintronlen G,-8.0,1.0 --min-intronlen 20 --max-intronlen 500000 --threads 12 -x rn6/genome -q -1 YH2mini_R1.fastq.gz -2 YH2mini_R2.fastq.gz |  samtools sort -m20G -l 0 -T /Users/me/Documents/TEMP/b -O bam | samtools view -O bam -@ 12 -o YH2mini.bam

Stringtie:

stringtie -fr -p 12 -G rat6.gtf -e -b ballgown/YH2mini -f 0.15 -m 200
-a 10 -j 1 -c 2 -g 50 -M 0.95 -o ballgown/YH2mini/YH2mini.gtf YH2mini.bam

PrepDE:

prepDE.py -i ballgown -g YH2mini_gene_counts.csv -t
YH2mini_transcript_counts.csv -l 100

So peep can see I'm not lazy, I have a work around where I can run this, extra prior to prepDE.py, but seems very wrong!!
if(CORRECTGTF){

        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/",PREFIX,".gtf > temp.txt")
        system(command)
        command <- paste0("sed 's/transcript\\://g' temp.txt > ballgown/",PREFIX,"/",PREFIX,".gtf")
        system(command)
        command <- paste0("sed 's/gene\\://g' ballgown/",PREFIX,"/","t_data.ctab > t_datafix.ctab.tmp")
        system(command)
        command <- paste0("sed 's/transcript\\://g ' t_datafix.ctab.tmp > ballgown/",PREFIX,"/","t_data.ctab")
        system(command)
    }

The gene:ENSRNOG00000046319 in the GTF file is pretty standard. I don't think we can help without access to the actual data. Can you give us some examples? for instance, you seem to be searching through something for the gene names. What are you searching through? Presumably the output of your pipeline, but what kind of output is that? Then, you say that you sometimes see both but you show a count of 0 for gene:ENSRNOG00000004937. My guess is that you are getting the gene names from different fields, but I would need to understand what you're doing in more detail to be sure.

@ Terdon Thanks, I will just need to investigate now how to share data and upload anything needed (after a meeting now!).

We don't need real data, per se, just some examples to understand what you mean.

...OK, sorry... don't really know how to give examples beyond the one above and so the actual code and data are here... https://www.dropbox.com/sh/lqk0usgt8n2suvx/AADqZDNpKKPKFGWZQ9sn9ix6a?dl=0


text:
How to identify sequence origin in DNA shuffle reads?
The data
I have a large number of reads from sequences that were generated by randomly shuffling regions of two parent sequences together. See the following image:

The regions that are shuffled are around 100-200 bp in length. The full shuffled sequence is around 1000 bp long. In general, the regions will be in the same order as the original sequence. The shuffled sequences are fragmented and sequenced into 150bp reads. One read can contain more than two borders between shuffle regions from a different parent sequence.
The sequences are very similar (95%). Because of this high similarity you cannot pinpoint the border between two regions. The border will be a range between a variant that is unique to source1 and a variant that is unique to source2. That is ok because I am not interested in knowing from which source an identical region came.
The question
I would like to identify which part of the read originates from parent sequence one or two.
What would be a suitable tool to do this?
My attempts
I tried to look at some multiple sequence alignment tools like MAFFT or Clustal but I am not interested in aligning all reads to each other. They also don't identify the shuffle regions well in my testing.

This is the sort of thing that LAST has been specifically designed to work well for:

LAST is designed for moderately large data (e.g. genomes, DNA reads,
proteomes). It's especially good at finding rearrangements and recombinations: we believe last-split does that more rigorously than anything else.

Here's my usual approach for mapping with LAST:
lastdb -uRY4 -R01 reference_sequences.fa reference_sequences.fa
# training is probably not needed if your sub-sequences have no
# variation from the parental sequences
last-train reference_sequences.fa input_sequences.fa > trained.mat
lastal -P 10 -p trained.mat --split reference_sequences.fa input_sequences.fa > aligned.maf

LAST is very customisable. If you want to adjust the mapping or split sensitivity, have a look at the lastal documentation here, and the last-split documentation here.
For example, there's a mismap probability adjustment which excludes alignments with a mismap probability above a certain threshold:
# only form alignments if the mismap probability is 0.1% or less
lastal --split-m 0.001 ...

I tried a lot of the parameters. The mismap threshold did not have an impact for me. I found that lowering the minimum alignment score with the parameter lastal -e to 50 worked well to increase the sensitivity. It is now working perfectly.


text:
Number of Spots in a Spatial Seurat object
I have got an integrated seurat object of 21 Spatial samples. I want to know how many spots are in all the samples together, Is there a way I can do that? Also, Please suggest any R Packages that do good ploting with spatial data (Like dittoSeq for Single cell).
Thanks a lot,
AMK


text:
Is it common to get different number of SNVs+Indels across samples from vcf files generated using GATK and DRAGEN (counts are higher for GATK)?
We have 2 vcf (Whole Exome Sequencing (WES) data; germline samples) files (e.g., vcf_1 and vcf_2). vcf_1 was generated (Ref. genome: hg38) using the GATK pipeline for 250 children and their parents and vcf_2 was generated (Ref. genome: hg38) using the DRAGEN pipeline for another 50 children and their parents (Samples are different). These are raw vcf files and I have not applied any quality-related false positive filters to them.
Now, when we compare the sample-wise SNV counts from vcf_1 and vcf_2, there is a huge difference. vcf_2 contains on average 4-6 times lesser number of SNVs per sample than vcf_1.
N.B:
I found this article (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4418901/), where it has been reported that the mean numbers of single-nucleotide variants (SNVs) and small insertions/deletions (indels) detected per sample was around 98k, respectively, for WES.  For our case, in the case of GATK-generated vcf samples, this varies from 130k to 188k, and for DRAGEN-generated vcf samples, this varies from 30k to 40k.
I understand both pipelines (GATK and DRAGEN) implement different methods to generate vcf files but is it common/normal to examine that much difference in SNV counts? I would appreciate it if someone can give some explanation (or provide some materials) about it.

We need a lot more detail to be able to help. What algorithm was used? What version of the various tools? Are these germline samples or somatic? What species? What kind of sample (WGS, WES, panels)? What filtering options were applied? Please [edit] your question and add more detail.

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.

I have updated the question with a few details (WES data, human sample, Ref. genome used: hg38, etc.) as you have mentioned. As soon as I will get some more info about the commands used to generate the vcf files or the version of GATK and DRAGEN, I will update the question further. Many thanks in advance.

I suspect what's happening here is an unfortunate side-effect of making variant-finding algorithms memory-efficient, so that they can work on extremely large populations. By doing the processor-intensive variant calling separately on each individual, calculations can be easily parallelised to run quickly.
For the numbers you have (i.e. 300 children and their parents), there's an alternative approach for generating variant calling that uses all mapped files as input, rather than only one at a time, i.e. instead of running:
bcftools mpileup [options] in1.bam
bcftools mpileup [options] in2.bam
bcftools mpileup [options] ...

The all-at-once approach has a single mpileup call:
bcftools mpileup [options] in1.bam in2.bam ...

This will call multiple samples at the same time, so that the output will have identical variant lines for all individuals. It also has the added benefit of getting more information on variation, allowing more rare variants to be identified.
I'm not sure which approach is being used for DRAGEN vs GATK, but I do know that at some stage in the past the GATK standard practice was calling one sample at a time. Curiously, if this were still the case, I'd expect the results to be the other way round, as the one-at-a-time approach tends to be less sensitive.
I guess it's also possible that the pipelines have different sensitivity / specificity tweaks, that could also explain the differences in variant counts.


text:
From RNA-seq raw counts to normalized data
I have several RNA-seq datasets. Some of them provide RNA-seq raw counts, some provide FPKM, RPKM and some have transcripts per million (TPM) data. Non of them provide fastq files, all data is processed already. At the end I want all datasets to be normalized to TPM.
I'm using this code in order to normalize raw counts to TPM: (using R)
rpkm <- apply(X = subset(dataset),
                MARGIN = 2,
                FUN = function(x) {
                  10^9 * x / genelength / sum(as.numeric(x))
                })

TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

When RPKM is provided, and no raw counts is available, I use the second line in the same code:
TPM <- apply(rpkm, 2, function(x) x / sum(as.numeric(x)) * 10^6) %>% as.data.frame()

And when FPKM is provided, I use this formula to transform the data to TPM:
TPM = FPKM*X

where X = 1e6/[sum of all FPKM of a sample]
NOTE: genelength is obtained using the biomart package in R, to get the transcript length directly out of ensemble.
However, because I know the steps of TPM normalization in theory, one should firstly normalize to gene length, and then to gene depth. I'm not sure this is what the code does, feeling very skiptical about it even thought this code was given to me by a Phd student.
Can you guys please help me with this question? is my code right or should I alter something?
Thank you!

Regarding the use of TPM, it's a good idea to be aware of the issues associated with it, even if it's your only available option. There are a few different TPM questions that have been asked previously, this answer has a reasonable summary.

However, because I know the steps of TPM normalization in theory, one
should firstly normalize to gene length, and then to gene depth.

No, you should normalize to transcript length.  Which requires you to know the relative abundance of all the transcripts so you can correct for all their lengths.
I do not think it's possible to un-normalize RPKM or TPM without knowing all the lengths involved in doing the calculations.
Let me give an example: here is RSEM output for just one gene.  Note that the software thinks there are three transcripts present for this one gene:
transcript_id  gene_id length  effective_length    expected_count  TPM FPKM    IsoPct
ENST00000373020 ENSG00000000003 3768    3718.54 179.95  13.88   7.94    48.61

ENST00000494424 ENSG00000000003 820 770.54  0.00    0.00    0.00    0.00

ENST00000496771 ENSG00000000003 1025    975.54  11.94   3.51    2.01    12.29

ENST00000612152 ENSG00000000003 3796    3746.54 0.00    0.00    0.00    0.00

ENST00000614008 ENSG00000000003 900 850.54  33.11   11.17   6.39    39.10

ENST00000373031 ENSG00000000005 1205    1155.54 0.00    0.00    0.00    0.00

ENST00000485971 ENSG00000000005 542 492.54  0.00    0.00    0.00    0.00

So, what would you say the effective length should be for the gene level calculations, and would you get that right without having this transcript-level breakdown?

I do have the transcript lengths, as I mentioned above I got them from ensemble.

Do you have transcript level counts?

No.. I assume this is something you get when you have BAM files? Anyway, I don't have that. I was able to get the transcript length which is necessary for the normalization. Do you think the code here provided by Michael Love is better?
https://support.bioconductor.org/p/91218/#91256

No, bam files do not magically give you transcript level counts.  Software like Kallisto and RSEM and Salmon do, and they do all the calculations for you.  Transcript length is not enough.  Genes have multiple transcripts of different sizes, and you can have a mix of isoforms in your samples.  The software I mentioned will magically handle all this for you.  I don't know that it's possible for you to do the magic in reverse.  I'm sure you can wave your hands and do some math and get numbers.  I don't know that those numbers will be accurate.

I see. According to Michael Love, as you said,  there is no way to calculate accurate TPM values out of count matrix. As I mentioned earlier I have no access to BAM files unforunatlely. But in this link I think I got my answer, they discuss this matter, and Michael Love provides here a simple code to normalize TPM out of count matrix. Thank you for your help and explenations. This is the link https://support.bioconductor.org/p/91218/#91256


text:
Reduce Overfit by removing insignificant proteins through PSEA analysis
I'm training a ML model for disease prediction using protein composition as input, but overfit is present. While looking to remove proteins and reduce multi-correlation, as remove composition predictors through PCA have worked but not enough, this problem appeared and managed to use Metascape (1) for it.
I'm using the Protein Set Enrichment Analysis (PSEA, GSEA for genes) for significance of function as is a method to identify classes of genes or proteins that are over-represented in a large set of genes or proteins, and may have an association with disease phenotypes (2). So:

A higher significance coefficient (log10pvalue) between two proteins suggest they have a similar function so one should be removed, or a lower coef shows the need to erase this quasi-independent protein from the analysis?

Is this analysis good enough to reduce overfitting by remove redundant observations?

References:
(1) https://metascape.org/gp/index.html#/main/step1
(2) https://www.sciencedirect.com/science/article/pii/S2376999819300406

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.

Overfitting is usually resolved by adding more data and not removing data sets. The technical blurb is 'augmenting' your training data set to revolve overfitting.
It is possible there is a very large multi-gene family involved distorting the classification, thus reducing that could reduce the overfit and cause the algorithm to spread the classification over a much broader group of genes. The general advice is to find an alternative genomic source of proteins and add that, not for this specific question but in reducing overfit per se.
There are technical issues such as:

a parameter split, does the training have this? If not its strongly recommended - depending on the algorithm.
algorithm, for example decision trees are known for overfitting*, whilst naive Bayes would find difficulty overfitting**.

I personally think multiple measures would resolve this. However, general advice is to look at alternative methods and training approaches before looking at slimming down a data set. I would suggest there needs to be good evidence - which would need to be formally presented. You could be 100% correct (probably are), but its not an orthodox approach.
Clearly there is information here that isn't presented. What you've done is a feature selection based on genes and we don't know what that result was nor the gene group that emerged from that analysis - so it is difficult to dissect this precisely. It could be there is a compelling biological reason for concern of overfitting based on the result, but we'll not know that. What I am certain about is the orthodoxy in resolving overfitting is not leveraged here and any data scientist would be aware of this.
*, this is why random forests are considered strong
**, there are a lot of different approaches to ML and how they are classifying is very different, some are more prone to overfitting than others. Naive Bayesian and 'linear regression' are the least sensitive to overfitting.

Please be aware I do ML within the Python framework and its one of the few areas I'm not aware of the ease of cross-over into other platforms. This is simply because its an analysis where Python excels.


text:
Correct way of merging scRNAseq datasets of healthy and tumor cell?
While merging datasets from similar biological conditions, but different experimental conditions is a well-studied topic, and there are many batch correction techniques available. What I want to understand is the right way to merge scRNA seq counts data of healthy and tumor cells.
Does individually log-normalising the counts, followed by a simple concatenation of the objects works, or there is a better way to do this?

I'm aware that you've tagged this as python, but unfortunately I'm most familiar with doing single cell analysis using R. Perhaps there's a python way to do something similar.
Seurat has recently included a workflow for reciprocal PCA integration, which the authors claim should work well when the combined datasets are more different than usual:

RPCA-based integration runs significantly faster, and also represents a more conservative approach where cells in different biological states are less likely to ‘align’ after integration. We therefore, recommend RPCA during integrative analysis where:

A substantial fraction of cells in one dataset have no matching type in the other
Datasets originate from the same platform (i.e. multiple lanes of 10x genomics)
There are a large number of datasets or cells to integrate (see INSERT LINK for more tips on integrating large datasets)

A tutorial workflow can be found here:
https://satijalab.org/seurat/articles/integration_rpca.html

Indeed, I am a python user, but this is already helpful. Thanks :)


text:
Clustering individuals by gene presence/absence
I have a binary matrix with individuals as row names, and gene name as column name.
So, if the gene is present in an individual, we have 1, otherwise 0.
I would like to cluster individuals based on the gene presence/absence pattern.
Is it right to perform an MCA? Or a logistic PCA is needed?
Thanks,
Marco

Would you consider hierarchical clustering, i.e. making a tree? However you need to present a test data set and I could give you the tree with an explanation. PCA is okay, but trees have advantages, for reasons I could explain.

I believe the answers to your question are here:
https://help.ezbiocloud.net/gene-presenceabsence-analysis/
https://www.biostars.org/p/66355/


text:
Get genome coordinates from UniProt entry?
I'm trying to get genome coordinate information programmatically for multiple genes.
As an example, I'd like to get the annotation (list of exons and their chromosomes, starts, and ends) for this UniProt entry: https://www.uniprot.org/uniprotkb/B6TU39/entry.
There is a corresponding GenBank entry, which gives the CDS (https://www.ncbi.nlm.nih.gov/nuccore/EU968504) but I can't see anywhere how to extract the genome coordinates. How can I do this?
Thanks!

You have asked for a 'programmatic' solution; what programming language?

I'll be using Python, but anything using the UniProt/GenBank/etc. REST APIs would work.

In the first instance it'd be useful just to know how to do it using the website GUIs, as that's likely to give enough clues

Talking of web GUIs, you can use Ensembl plants. Put the uniprot accession in the search box and you will find the gene page.

You don't want the EU968504, that's the mRNA not the genome sequence. You want https://www.ncbi.nlm.nih.gov/nuccore/CM007647, also linked to from the UniProt page, just under the link to EU968504:

You also need the gene name, ZEAMMB73_Zm00001d029604 in this case, available further up on the page in the "Names & Taxonomy" section:

Next, we can download the GenBank file for the relevant chromosomal contig  using NCBI's entrez API:
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=CM007647&rettype=gbwithparts&retmode=text" > CM007647

This will take a few minutes. Once that is done, you can parse the file and extract the start and end coordinates using something like this:
#!/usr/bin/env python3
from Bio import SeqIO
import sys

gb_record=SeqIO.read(open(sys.argv[1],"r"), "genbank")

for feature in gb_record.features:
    if feature.type == "gene" and sys.argv[2] in feature.qualifiers['locus_tag']:
        if feature.location.strand == 1:
            strand = "+"
        elif feature.location.strand == -1:
            strand = "-"
        else:
            print("Unknown strand: ", feature.location.strand)
            sys.exit(1)

        print("%s\t%s\t%s\t%s\t%s\t%s\t%s\t.\t." % (gb_record.id, "GenBank",
                                                    feature.type, feature.location.start,
                                                    feature.location.end, ".", strand))

You pass the GenBank file and gene name as arguments, and this is the result:
$ foo.py CM007647  ZEAMMB73_Zm00001d029604
CM007647.1  GenBank gene    78770080    78773320    .   +   .   .

Thanks so much for the detailed answer!

@MichaelDunne CM007647 comes from GRAMENE-4.0, not the latest NAM-5.0 maize genome. On NAM-5.0, this gene is located at 1:77,982,020-77,986,746.

You could align back the peptide sequence to the genome, respecting splicing, for example using miniprot from Heng Li.

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.

Hi and welcome to the site! We like answers to be a bit more comprehensive here. Could you expand your answer to explain how to use miniprot to do what the OP asks for? Also, you would need to explain how to go from a protein accession to its fasta sequence to then use that to align to the genome, as well as how to figure out which genome to align to and how to differentiate between homologs.


text:
Retrieval protein definition from accession ID
I have a list of protein accessions:
XP_025107904.1
PVD22034.1
AJA37836.1

How can I obtain the definition of the protein, and the related organism?
Output:
XP_025107904.1  receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata

Thanks for the help,
Marco

I solved with a Python script (through biopython):
from Bio import Entrez
from Bio import SeqIO
handle = Entrez.efetch(db="protein", id= input_file, rettype="gb", retmode="text")
for seq_record in SeqIO.parse(handle, "genbank"):
print(seq_record.description + "\t" +  seq_record.id)

great! Could you please post your script as an answer and accept it so the question can be marked as answered?

@Marco could you please consider upvoting/accepting the answer. It does help with site stats and its only fair to terdon's answer

You can do this using the Entrez API. For example, using the three identifiers from your example, you can do:
curl "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=protein&id=XP_025107904.1,PVD22034.1,AJA37836.1" > file.xml

This will bring an XML file with the following contents:
$ cat file.xml
<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE eSummaryResult PUBLIC "-//NLM//DTD esummary v1 20041029//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20041029/esummary-v1.dtd">
<eSummaryResult>
<DocSum>
    <Id>1397721498</Id>
    <Item Name="Caption" Type="String">XP_025107904</Item>
    <Item Name="Title" Type="String">receptor-type tyrosine-protein phosphatase kappa-like isoform X2 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1397721498|ref|XP_025107904.1|[1397721498]</Item>
    <Item Name="Gi" Type="Integer">1397721498</Item>
    <Item Name="CreateDate" Type="String">2018/06/07</Item>
    <Item Name="UpdateDate" Type="String">2018/06/07</Item>
    <Item Name="Flags" Type="Integer">512</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">1271</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">XP_025107904.1</Item>
</DocSum>

<DocSum>
    <Id>1383748625</Id>
    <Item Name="Caption" Type="String">PVD22034</Item>
    <Item Name="Title" Type="String">hypothetical protein C0Q70_17837 [Pomacea canaliculata]</Item>
    <Item Name="Extra" Type="String">gi|1383748625|gb|PVD22034.1||gnl|WGS:PZQS|C0Q70_17837[1383748625]</Item>
    <Item Name="Gi" Type="Integer">1383748625</Item>
    <Item Name="CreateDate" Type="String">2018/04/26</Item>
    <Item Name="UpdateDate" Type="String">2018/04/26</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">400727</Item>
    <Item Name="Length" Type="Integer">602</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">PVD22034.1</Item>
</DocSum>

<DocSum>
    <Id>735659541</Id>
    <Item Name="Caption" Type="String">AJA37836</Item>
    <Item Name="Title" Type="String">fucolectin-related molecule [Littorina littorea]</Item>
    <Item Name="Extra" Type="String">gi|735659541|gb|AJA37836.1|[735659541]</Item>
    <Item Name="Gi" Type="Integer">735659541</Item>
    <Item Name="CreateDate" Type="String">2014/12/25</Item>
    <Item Name="UpdateDate" Type="String">2014/12/25</Item>
    <Item Name="Flags" Type="Integer">0</Item>
    <Item Name="TaxId" Type="Integer">31216</Item>
    <Item Name="Length" Type="Integer">722</Item>
    <Item Name="Status" Type="String">live</Item>
    <Item Name="ReplacedBy" Type="String"></Item>
    <Item Name="Comment" Type="String"><![CDATA[  ]]></Item>
    <Item Name="AccessionVersion" Type="String">AJA37836.1</Item>
</DocSum>

</eSummaryResult>

You can now parse this in your favorite language/parser and get the information you need. I... am not too good with XML so my attempts with xmllint failed, but here's a naive awk approach:
$ awk -F'[][><]' '{ if(/Name="Title"/){desc=$3; spc=$4}else if(/Name="AccessionVersion"/){name=$3} if($0=="</DocSum>"){print name,desc,spc; name=""; desc=""; spc=""}}' file.xml
XP_025107904.1 receptor-type tyrosine-protein phosphatase kappa-like isoform X2  Pomacea canaliculata
PVD22034.1 hypothetical protein C0Q70_17837  Pomacea canaliculata
AJA37836.1 fucolectin-related molecule  Littorina littorea

I stress that this is a very bad idea and you should use a proper XML parser, but it might work...


text:
subset from a dataframe using other dataframe in R
I have a data frame (named as df) with single column of protein IDs and other file with two columns of protein IDs (named as df1). I need to find the df IDs present in either first or second column of df1 and extract corresponding Id from the df1.
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612","ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000160293","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000142949","ENSG00000112936","ENSG00000162434"))

Sample output:
df2= data.frame(p = c("ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400","ENSG00000198400"),p2 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119","ENSG00000112936","ENSG00000162434"))

This seems to be an XY problem. Can you give us an idea of your bigger picture problem?

I'd say this belongs on Stack Overflow as is more of a pure coding problem - and you'll get a faster response there. I'd suggest deleting this and reposting here.

Gene names... looks like there's a biology aspect to this; adding more context to the question would make it more obviously bioinformatics.

you can try this:
df = data.frame(p = c("ENSG00000198400","ENSG00000066044","ENSG00000188612",
                       "ENSG00000177606","ENSG00000146648"))
df1 = data.frame(p1 = c("ENSG00000105647","ENSG00000110395","ENSG00000165119",
                         "ENSG00000160293","ENSG00000198400","ENSG00000198400"),
                 p2 = c("ENSG00000198400","ENSG00000198400","ENSG00000198400",
                         "ENSG00000142949","ENSG00000112936","ENSG00000162434"))

df2 = rbind(df1[df1$p1 %in% df$p,],df1[df1$p2 %in% df$p,])

df2 <- transform(df2, p2 = ifelse(df2$p2%in%df$p, p1, p2), p1 = ifelse(df2$p2%in%df$p, p2, p1))


text:
Comparison of fastq files reads
My goal is to compare reads from two different fastq files on a Linux machine.
The following are the comparisons to perform:

How many common reads are between the two fastq files?
How many reads are present in one fastq file that are not present in the other? i.e reads that are unique to each file.
Where and how do the sequences differ in both files? Is there a difference in nucleotide base at the end of the sequence(3' end), the start of the sequence(5'end) or is it somewhere in the middle of the sequence? Is there a difference in length of the sequence? If so, what is the length and what is the difference? I would like to get the original sequences plus the difference between them for analysis.
Does the quality score differ? Where and how do the quality of sequences differs?

Any guidance regarding how to approach this question is very appreciated.
Edit: I'm working on paired-end reads from the same individual. What I want to do is to compare different trimming tools from their output.

Are you expecting the files to be similar? Are these WGS samples? Targeted? Using the same assay? Are they reads from the same individual? Same species but different individuals? Different species? The details will depend on how many differences you expect. Also, this is a whole project, what do you have so far, and what do you need from us?

Are you absolutely sure that comparing the files at this level is the right way to go about answering whatever questions you have?  Fastqs have millions and millions of reads, you really want to find the differences on a letter by letter basis?

Please show example input files (with enough lines to make the problem obvious), and example output

https://github.com/OpenGene/fastp

I doubt anyone has time nor energy to write full code for you to do this but I would start with Bio.SeqIO from BioPython, you could start here: https://www.biostars.org/p/147317/

There are several things you can try. For example, Simka (https://github.com/GATB/simka) is tool that compares the profiles of the genomic content of fastq files and generates several types of plots e.g. heatmap.
Another tool that might be useful is multifastqc (https://multiqc.info/), which aggregates results and has multiple options.


text:
Snakemake MissingRuleException
I'm trying to run a snakefile for the first time with limited coding experience using salmon to index a reference genome. I'm not too sure what I'm doing wrong based on this error message.
Error: Building DAG of jobs...
MissingRuleException:
No rule to produce snakefile (if you use input functions make sure that they don't raise unexpected exceptions).

Code:
rule salmon_index:
        input:
                R1="sequences/*R1.fastq.gz",
                R2="sequences/*R2.fastq.gz"
        output:
                Index="results/salmon_index"
        shell:
                """
                salmon index --threads 8 -t references/genome.fa.gz -i salmon_index
                """

You're not using wildcards the way snakemake is equipped to work with them. See:

FAQs such as https://snakemake.readthedocs.io/en/stable/project_info/faq.html#how-do-i-run-my-rule-on-all-files-of-a-certain-directory and https://snakemake.readthedocs.io/en/stable/project_info/faq.html#i-don-t-want-expand-to-use-the-product-of-every-wildcard-what-can-i-do
https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#wildcards