text:
Methods to work directly with bigbed file?
I am now trying to find intersect genes of two files: one is .bigbed file (very large, 117GB), and the other is a .bed file.
At first, I tried to convert bigbed file to bed file and then applied bedtools. I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server. Thus I planned to convert the .bb file on my own computer and then upload the converted .bed file to the server, and then apply bedtools. Firstly I tried to convert some smaller .bb data sets (< 10 GB) to .bed and it worked slowly but well on my computer. However, when I tried the 117GB data there was an error:
Invalid argument, Error reading 9188986873 bytes.
My code is :
bigBedToBed JASPAR2022_hg19.bb hg1 9.bed
Thank you in advance!

Now the bug is solved, and I got many other good suggestions. Below you can find many methods to work with the bigbed file, directly working with it or converting it to bed file.

You asked

Does anyone know any tools that can intersect two bigbed files

...but it appears you really wanted to subset a single very large bigbed file based on the contents of a much smaller bed file.
I disagree with your summation "there is not a very efficient way to work directly with the bigbed file".  I see two options.
The first is to repeatedly call bigBedTobed once for each row in your .bed file.  If you have multiple cores, this can be parallelized using GNU parallel.  The following does that, emitting each overlapping result with additional column being the index into your .bed file from which the result stems.
cat your.bed | parallel --lb -j --colsep '\t' bigBedToBed -chrom={1} -start={2} -end={3}  your.bb  /dev/stdout '|' sed -e 's/$/\\t${PARALLEL_SEQ}/'

Depending upon your application, there is a second approach which uses R/BioConductor's rtracklayer::import.bb which can be parallelized across (chunks of) your .bed file.
I considered these approaches in what I expect is an almost identical situation:
My bigbed holds genome-wide transcription factor binding site predictions from JASPAR Genome Browser tracks as produced by wassermanlab / JASPAR-UCSC-tracks.
My bed file holds the loci of ~23K putative enhancers as determined by an analysis of regions of accessible chromatin from ATAC-Seq and various ChIP marks.
In my case I wished to produce a sparse Matrix with one row for each putative enhancer, and one column for each JASPAR motif (which will form the basis of downstream enrichment analyses).
bbf<-BigBedFile('path/to/your.bb')
gr<-import.bed('path/to/your.bed')
## I found I had to harmonize seqlevels (e.g. chromosome identifiers) to make this work.  YMMV
gr<-renameSeqlevels(gr,paste0('chr',gsub('^(K.*)\\.(\\d+)$','Un_\\1v\\2',seqlevels(gr))))
mm<-motifMatrix(gr,bbf)

which results in a sparse Matrix, viz:
> rownames(mm)<-as.character(gr) # optional.
> mm[1:5,1:5]
5 x 5 sparse Matrix of class "dgCMatrix"
                 MA0002.2 MA0003.4 MA0004.1 MA0006.1 MA0007.3
chr1:11156-12180       14        .        2       28        2
chr1:17836-20023       25        .       27       63        2
chr1:22495-24367       22        1       22      137        2
chr1:27533-28038        2        2        4        7        .
chr1:36549-37079        6        .        2       41        .

This depends upon the following definition:
motifMatrix<-function(gr,bbf,motifID=c('name','TFName')[[1]],mc.jobs=100) {
  ## produce a sparse Matrix with one row for each region in
  ## GenomicRanges gr, and one column for each JASPAR motif indicated
  ## as overlapping <gr> accorinding to BigBedFile <bbf>.  Default to
  ## using the 'name' column as the motif identifier, but allow for
  ## possibly using another, such as 'TFName' (as used by 2022 version
  ## of
  ## https://jaspar.genereg.net/genome-tracks/#ucsc_tracks). Parallelize
  ## into mc.jobs queries, defaulting to 100, which was found suitable
  ## for ~25000 loci in <gr> against JASPAR matches to zebrafish
  ## danRer11 genome running on a few dozen cores.  In general,
  ## knowing how "big" each parallel job *should* be depends on
  ## available RAM and cores, and number of loci in <gr>.
  ##
  ## AUTHOR: malcolm_cook@stowers.org
  ##
  ## NB: you must install a version of `rtracklayer` which addresses
  ## issue [import\.bb returning incorrect ranges & results as
  ## compared with kent tools
  ## bigBedToBed](https://github.com/lawremi/rtracklayer/issues/59)
  ## (such as by `install_github('lawremi/rtracklayer')`)
  si<-parallel::splitIndices(length(gr),mc.jobs)
  mm<-mclapply(si,function(i) {
    gri<-gr[i]
    m<-import.bb(bbf,selection=BigBedSelection(gri, colnames = c(motifID)))
    fo<-findOverlaps(gri,m)
    xt<-
      xtabs(~.,data=list(locus=i[queryHits(fo)],motif=mcols(m[subjectHits(fo)])[[motifID]]),sparse=TRUE)
    xt
  })
  mm<-
    ## rBind.fill(mm) ## NOT!  It does NOT preserve dimnames and
    ## returns matrix when given Matrix.  Issue reported:
    ## [[https://github.com/cvarrichio/Matrix.utils/issues/5][rBind.fill
    ## on a list of Matrix arguably should return Matrix but returns
    ## matrix #5]]
    Reduce(rBind.fill,mm) # which does preserve dimnames and return Matrix
  mm
}

NB: If you choose this approach, you MUST install a version of rtracklayer which addresses issue the import.bb returning incorrect ranges & results as compared with kent tools bigBedToBed (such as by install_github('lawremi/rtracklayer'))

Yes Malcook, this is exactly what I want at first! Thank you so much!

glad to help - would you consider to change the accepted solution to mine?

This is the accepted answer for sure @BinhuanSun

No matter how you approach this problem 110+ GB is big:
Your options are:

Break up the file and perform the intersection on a chunk by chunk basis
Ask your system administrator to install the program
Load it into Docker, your systems administrator will understand this. Installing docker is a single line of code (which I forget). Docker is a container which protects the system from the program and no systems administrator in this day and age would refuse a request for Docker installation.

Your.systrm administrator might be concerned about the file sizes involved, but I think a better dialogue with them is the right way forward here.

the first option is maybe a good one, you can select individual chromosomes at a time with bigBedToBed

Thank you for your suggestions! I will try the first one on my computer first :)

@ColinD Yes, I will try. Thank you :)

Hello M, because today is Sunday so I first tried the first method you suggested. I tried to break up the file using argument -Chrom of the bigBedtobed program: bigBedToBed -chrom=chr1 JASPAR2022_hg19.bb test.bed. However, I still got the same errorInvalid argument, Error reading 9188986873 bytes. I then tried a relatively smaller data (66GB) but got the same error. It seems that on my computer the bigBedtobed program only works on data that is smaller than 10 GB. I still have over 600 GB of drive memory, so I think this may due to the lack of RAM?

If its been written in C (I honestly dunno) a core-memory dump is a common error. However, if you > 10 chromosomes each file should be smaller than 10 GB? If its Drosophila then you get the same error. Could you check the file size for each intersect?

Hello M, Thank you for your quick reply! Actually, it is a Genomic track and TFBS prediction for human genome, and it is a whole file. I think after breaking up every file will < 10GB, but the thing is I don't know how to split the .bb file directly and then convert them to the .bed file :(. The website for the file I am working with is here: link. The file name is JASPAR2022_hg19.bb.

I think this should be a separate question. It wouldn't be hard to modify the file to allow breakup the file given the chromosome tag is a recognised tag for this purpose as @ColinD has stated. Perl pie will modify massive files with ease (if it was simplistic).

Yes M. I think you are right, and today I asked my supervisor to help me apply for an account on another server that installed those UCSC programs, and I hope it can work on the server. Thank you again for the patient reply :)

The question I asked a couple of days ago has been solved. Here I want to make a summary of this and hope it can help people who also encounter such problems later.
First, about the way to process the bigbed file. I didn't find a very good way to work directly with bigbed file at that time, so I still firstly converted the bigbed file to the bed file using bigBedTobed program and then apply bedtools. The bigBedTobed program can be downloaded through this link: UCSC programs, and it also contains many useful programs from  UCSC genome browser. I finally converted my data (117 GB bigbed file) on the server, and it took ~2.5h to convert it to bed file. After conversion, the converted bed file is about 550 GB. And it took ~4h to intersect with a 4 Mb bed file using bedtools.
Second, the bug that I encountered when I tried to convert the bigbed file on my own laptop. The code I used is bigBedToBed JASPAR2022_hg19.bb hg1 9.bed and it has error Invalid argument, Error reading 9188986873 bytes. I reported this error to the UCSC genome browser group and they said it is a MacOS specific bug, reaching a default data limit. They have solved it in their development environment but it still needs some time to reach the public site. The current solution for this is to use a URL to the file instead of a local file allows the program to read appropriately sized pieces of data at a time(-udcDir= in bigBedToBed).
Finally, this is my first time asking questions here and I really appreciate those good suggestions from @M__ and @ColinD.

Renew: Methods from @malcook can work directly with the bigbed file! For a very large bigbed file, it is indeed a more efficient way.

Hi nice to hear from you @BinhuanSun. Cool glad its worked out and thanks for the feedback.

MacOSX bug .. thanks thats interesting, I encountered something similar. They are nice machines, but it ain't a cool way to do things IMO.

My 2p, a little late and maybe not too useful:

I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server.

Presumably, you cannot install in /usr/local/bin or some other root directory because you don't have admin rights. This is the common setup on shared servers. However, many programs can be installed in the user's home directory without admin rights. Even better, install conda and mamba first (no need to be admin) and use that to install other programs.
The problem with docker is that some system administrators are concerned about security and refuse to install it (and docker does need admin rights). Whether such concerns are founded, I don't know...

I tried to convert bigbed file to bed file and then applied bedtools, but I found it was very slow and couldn't make it.

bedtools intersect supports the -sorted options that should speed things up quite a bit - Have you applied it? Your huge bigbed files should be already sorted after you convert to bed.

Hi Dariober, yes later I found I could install programs in my home directory. Sorry for this super easy problem because I am not familiar with the Linux server. I just checked my converted bed file, and it is indeed sorted, so I think it is just because the converted bed file is too big (560GB).


text:
Chimera- How to go about changing the threshold of the EM map?
It is my first time using Chimera. I have uploaded a PDB file on to Chimera and opened it up on the software. I need to change the threshold of the EM map so that I can analyze what the EM map would look like with a threshold too high and too low.
Please let me know if you have any suggestions.
Thanks

Could you supply more details on the goal of your analysis?

https://www.cgl.ucsf.edu/chimera/current/docs/ContributedSoftware/volumeviewer/framevolumeviewer.html

I am trying to view a protein with a threshold that Is too high, a threshold that is too low, and finally, the correct threshold to se the secondary structures.

source: Stack Exchange

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

sort: votes


text:
How to visualize electron density of a specific part of the protein residue?
Hi everyone (mostly who are interested in structural bioinformatics).
I face the problem that I cannot plot the electron density of a specific part of my protein. For the whole protein it is not a problem:
from Bio.PDB import *
import nglview as nv
import ipywidgets

pdbl = PDBList()
path = pdbl.retrieve_pdb_file('1FAT')
parser = MMCIFParser()
struct = parser.get_structure("1fat", path)
view = nv.show_biopython(struct)
view

But when I tried to pick up residues from 10 till 109 (res_list = Selection.unfold_entities(struct, "R")[10:109]) and plot I can't do that, because for nglview I need to have a special Object. As I'm a newbie, I'd like to address this question to specialists.
At the moment I'm trying to visualise model (structure), not electron density. But the next step is trying to find a package to visualize electron density also. I've just started becoming familiar with some open packages.
Thank You in Advance!

Hi @KitsuneBreeze its a good question and looks to be resolved. Please remember to "accept" and "upvote"

NGLview runs off the JavaScript libray NGL, which is documented as a series of very helpful examples: in the gallery. The commands are not one-to-one (slightly frustratingly), but it is helpful as you can figure it out easily with a few dir and help commands. In this case whatever example of density being loaded it will reveals that it's a different component —this is the same everywhere as density files do not contain atomic coordinates. Unfortunately, there's not a direct way to generate a new ComponentViewer from different sources to an NGLWidget and you cannot merge two of them: the various .show_* functions return an NGLWidget but aren't part of it.
First, I am not sure why you are going through a convoluted route for loading a PDB code:
pdb_code = '1fat'

import nglview as nv

view: nv.NGLWidget = nv.show_pdbid(pdb_code)

Now, on loading the component a representation is added to it.
These can hide the representation you then add and things will not make sense, until this detail is recalled! I fall in that pit nearly every time.
But we can change it with update command:
# I don't like wonky sheets
view.update_cartoon(smoothSheets=True, color='salmon', component=0)

Now let's add the density map the circuitous way (no .show_url).
import requests, io
response: requests.Response = requests.get('https://www.ebi.ac.uk/pdbe/coordinates/files/1fat.ccp4')
response.raise_for_status()
fh = io.BytesIO(response.content)
comp: nv.component.ComponentViewer = view.add_component(fh, ext='ccp4')

Now we can add a surface representation, but first let's clear the default.
For the arguments I had a gander at the example:
view.clear_representations(component=1)
view.add_representation('surface',
                        color='grey',
                        wireframe=True,
                        isolevel=1.5,
                        component=1)
view.center('100-110:A')
view

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


text:
When can we have duplicated variants in a VCF file, excluding artifacts?
I read these two questions on duplicated VCF variants: 1, 2.
However, it is not entirely clear to me when duplicated variants can arise (excluding merging artifacts). I want to understand if duplicated variants should be always considered artifacts, or if there are specific use cases where they are needed.

Do you have a specific problem in mind? There are multiple reasons why duplicate variants may appear (especially when merging datasets), and Stack Exchange doesn't really work well when there are multiple "correct" answers for the same question.

I don't have a specific problem in mind. I wanted to understand when to expect duplicate variants, since I noticed that some algorithms filter them. Although I understand merging artifacts, I wonder if there are additional corner cases where these duplicate variants appear. The accepted answer can be a community wiki where multiple "correct" answers can be collected, I believe.

Okay, probably better asking this question on Reddit r/bioinformatics then.


text:
Displaying BIC and AICs from MSA nucleotide multisequence fasta file
I am initiating a project reproducing MEGA style trees using R to find the best DNA protein models under AIC/BIC. I have read the R AICs and BIC documentation but I'm struggling to implement this.
Could I implement AICs and BICs without having to specify the number of sequences in the fasta file?
My initial code that works the same as MEGA for aligning and building a phylogenetic tree where I want to view AICs and BIC (as per reference picture):
library(seqinr)
library(adegenet)
library(ape)
library(ggtree)
library(DECIPHER)
library(Biostrings)
library(viridis)
library(ggplot2)
library(msa)

mySequenceFile <- system.file ("sequences", "allSequences.fasta", package ="msa")
mySequence <- readDNAStringSet(mySequenceFile)

mySequence

myAlignment <- msa(mySequence, "Muscle")

myAlignment

alignment <- msaConvert(myAlignment, "seqinr::alignment")
distMatrix <- dist.alignment(alignnment, "similarity")
clustering <- hclust(distMatrix)

plot(clustering)

dendrogram <- as.dendrogram(clustering)
phylotree <- as.phylo(clustering)

plot(phylotree, type="radial")

Please ignore some libraries as I previously ran into some "-" deletion issues when I was writing the code

Having looked carefully at the question, I understand the 'gap in the market' but all this code will do is generate a provisional tree. I would suggest looking at protein models first because you need models to perform an AIC/BIC. So far you haven't generated used or described a protein model and without that there's nothing to compare via AIC.
To answer the direct question you should seek to avoid removing the sample size from the AIC and base the work on the equation below. If you have large samples sizes n can be removed (which is what you are wanting).
The work you should read is jModeltest which is here, https://github.com/ddarriba/jmodeltest2.
This originated from Modeltest by David Pasado.
Their code on GitHub will tell you exactly how to implement AIC and BIC. I recommend ditching R in this instance. If done via subprocessing (below) its not hard.
Their work concerns nucleotides and never really hit the mainstream, but was used by some.
The key criticism is like your approach that the tree generated was based on neighbor-joining. The tree will influence the matrix so there is a danger that the ML or Bayesian phylogeny resulting from the model will be biased by neighbour-joining. Secondly, it was not computationally expensive to do this entirely by ML or Bayesian approaches because it's a 4*4 matrix.
Protein phylogeny is very different, obviously the matrix is bigger 20*20 the other issue is that it's a 'static matrix'. There are approaches where they "combine" matrices but I don't know the algorithm. Thus we know exactly how the protein matrices are generated, "combined" needs a clearer explanation and how this is obtained could be a good application for AIC. Thus could provide useful clarity. There is in addition.the potential to generate new matrices and assess them, but the concept of introducing meaningful dynamism into the matrix does have potential and this is absent (except for cryptic combined options). AIC would then assess whether this is helpful. It might not be hugely helpful in traditional static matrices, but would compare them against further hypotheses of mutation and innovative approaches. If used deployed correctly could be cool because it allows you to experiment and ask the question what precisely does this protein matrix mean for my data?
More generally stretching the utility of AIC type metrics is very interesting in 'bleeding edge' machine learning (admittedly its a bit futuristic).

I suspect (could be wrong) jModelTest is coded in Java, an alternative would be to subprocess via Python using relevant compiled code. This was how the original Modeltest was written, essentially driven through the leading maximum likelihood package of the time known as PAUP. Python does implement AIC/BIC via the statsmodel library, scikit-learn (machine learning) or use the equation directly.

k= number of independent variables to build model
L= maximum likelihood estimate of model
.. for small sample sizes (n)
R is not known for subprocessing and I am not aware of any ML/Bayes libraries for phylogeny.
In summary, this is a good project and it is worth investing the time because it has a clear, meaningful output and potential to extend into more innovative experimental approaches. I for one would use established  code (i.e. written by someone else) that implemented AIC/BIC for proteins matrices and would be fairly straight forward (for you) to write. How you maximise its use in deploying it to help understand protein matrices needs creative thinking (remember its parameters verse likelihood).
Good luck.

Maybe instead of using ML I can try to play around DL as well with some good datasets. Thanks for the the guidance, some things make more sense now and I will try to do everything using python instead.

Cool, DL ... deep learning? Scikit-learn does machine learning, not quite the same. ML/ AIC looks a good staging post (a reviewer will ask for BIC I'm afraid). Please do remember to upvote.

Yes I meant deep learning. I played around some python scripts last night and R as well. Will be creating a GitHub repository later on this week with my progress. Unfortunately I can't upvote yet due to my reputation but I choose it as the answer. Thanks again.

Deep learning isn't trivial and you'll likely encounter skepticism, ML sci-kit lean is widely accepted and much easier.


text:
Normalization methods for single cell RNA sequencing that take read count into account
I have two RNA-seq datasets.  One was sequenced at an average read count of 1.5 million per cell the other at 43K average reads per cell.  For the first I also have the meta data from reads alligned and unaliagnable for each cell, but I may have to do some digging for the other dataset.  I am looking for the presence of a few very informative genes being expressed at similar levels in the second dataset to confirmed examples of these cells in the first.  The first dataset is much smaller.
My question is this: Are there normalization methods which I could use given the data I have to allow me to accurately say these cells are in the 90th or 80th percentile of a a particular gene being expressed?

More reads/UMIs will lead to more counts in the matrixes, so you can normalize for sequencing depth without using the fastq/bams directly. Indeed this is common with metrics like CPM. However, what are the single cell technologies used in this case?

The default normalisation methods of DESeq2 have a quantile-like normalisation process (based on raw read counts), and this default process should be able to handle datasets with different read depths.
If that's not good enough, DESeq2 allows you to provide alternative information (e.g. cell counts) to bypass the usual normalisation process using the sizeFactors function:

The sizeFactors vector assigns to each column of the count matrix a value, the size factor, such that count values in the columns can be brought to a common scale by dividing by the corresponding size factor (as performed by counts(dds, normalized=TRUE)). See DESeq for a description of the use of size factors. If gene-specific normalization is desired for each sample, use normalizationFactors.

More details here:
https://www.rdocumentation.org/packages/DESeq2/versions/1.12.3/topics/sizeFactors


text:
Where can we download a curated dataframe of the full archive of UniProt?
I want to download a curated dataframe of the full archive of UniProt.
Any format as long as it is a table (.csv, .tsv, SQL etc.).
Solution 1: Download the .fasta and do the formating myself -> very slow + information incomplete.
Solution 2: Use API like BioServices -> still very slow but more information.
Do we have an already-curated dataframe containing as complete as possible the entry + metadata of UniProt that can be downloaded right away?

Answer from @matteo-ferla, converted from comment:
There are many tables in the downloads page such tables enumerating each protein in a species proteome with basic stats, but not all the data. The uniprot records are very structured so would convert poorly to a table (even filling cells with lists/arrays and dictionaries/objects)… so a fuller table would be a mess.


text:
Generating contrast matrix for limma in loop
I am trying to generate contrast matrix for different conditions using for loop.
cluster_len <- 3
req_contrasts <- c()

for (m in 1:cluster_len) {
  for (n in 1:cluster_len) {
    if (n > m) {
      req_contrasts <- c(req_contrasts, (paste0(paste0("Group", n), "-", paste0("Group", m))))

    }
  }
}

req_contrasts returns
"Group2-Group1" "Group3-Group1" "Group3-Group2"
The contrast that needs to be passed is
"Group2-Group1","Group3-Group1","Group3-Group2"
I have tried the following code to combine each string and form a single vector
# To separate each condition by a double quote and comma
temp_1 <- paste0(req_contrasts, collapse='","', sep="")
temp_1
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"
# To remove the escape character
temp_2 <- gsub("\\", "", temp_1, fixed = TRUE)
temp_2
[1] "Group2-Group1\",\"Group3-Group1\",\"Group3-Group2"

This does not work
contrast.matrix <- makeContrasts(temp_2, levels=design_mat)
Error in parse(text = ej) : <text>:1:14: unexpected string constant
1: Group2-Group1","
^

Limma throws the following error in presence of escape character. I am assuming it is because it considers that the string is yet to end.
Following input works
contrast.matrix_2 <- makeContrasts("Group2-Group1", "Group3-Group1", "Group3-Group2",                 levels=design_mat)

In my analysis the number of clusters to be analysed is 16. Hence, the contrast conditions cannot be generated manually.

I wrote this snippet a while back to generate all pairwise contrasts:

make_all_contrasts <- function (group, delim="_vs_"){

  suppressMessages(require(limma))

  #/ ensure that group levels are unique
  group <- sort(unique(as.character(group)))

  #/ make all combinations
  cb   <- combn(group, 2, FUN = function(x){paste0(x[1], "-", x[2])})

  #/ make contrasts
  contrasts<- limma::makeContrasts(contrasts=cb, levels=group)
  colnames(contrasts) <- gsub("-", delim, colnames(contrasts))

  return(contrasts)
}

Assuming three contrasts this will create a limma-like contrast matrix with all possible pairs:
make_all_contrasts(group=c("groupA", "groupB", "groupC"), delim="_vs_")

        Contrasts
Levels   groupA_vs_groupB groupA_vs_groupC groupB_vs_groupC
  groupA                1                1                0
  groupB               -1                0                1
  groupC                0               -1               -1

You can provide a function on combn to make the code easier to read: combn(group, 2, FUN = function(x){paste0(x[1], "-", x[2])}) and skip an unlist and a lapply call.

Cool, I did not know that, thanks!


text:
Comparing differential expression across samples - is batch effect correction needed?
I have a bulk RNA-seq dataset made up of control and treatment conditions for a range of cell lines. This dataset was generated in two batches, such that the cell lines are split between batches but both the treatment and control for each cell line are within the same batch. As the cell lines are very different, I'm looking to do DE analysis for each cell line and then compare these DE responses between cell lines.
I'm wondering whether batch effect correction is useful here since I'm only ever comparing DE values and not the raw counts themselves. Also, the baseline expression values for the cells lines are so different that they already constitute a strong biological batch effect that is perfectly correlated with the technical batch effect, which batch correction wouldn't be able to account for.

There is no need for batch correction because all your comparisons are within batch. As long as you are making comparisons for each line, and treatment and control for the same line are in the same batch, then any batch effects are already accounted. for.

If your treatment and control are prepared in different batches, it's impossible to distinguish between treatment-specific effects, and batch-specific effects. DESeq2 will complain when you attempt to carry out such correction in the statistical model for this reason.
Baseline differences in expression are a large reason for carrying out batch correction in the first place; it makes differential expression analyses more powerful, not less. I'd recommend consulting with a statistician (prior to carrying out an experiment / sequencing) regarding the best way to split samples into different batches so that batch effects can be discovered and accounted for. As you mention, there are multiple effects that should be considered (e.g. cell lines, sequencing run).

Could you expand a bit on why batch correction helps with baseline differences? For example, if I have multiple cell lines and treatments from a single experiment, sequenced in the same run, are you saying that I should still do batch correction to adjust for the baseline expression? Doesn't that run the risk of altering real biological differences?


text:
Is this a complete bacterial genome? (CheckM)
I've been working with a bacterial genome assembly. The initial stats look pretty good (3 main pieces, depth ~50x [20x - 100x] if I calculated it right). I'd expect this to be a virtually complete genome apart from some low-complexity regions.
However, when I use CheckM to estimate completeness / redundancy, I get ~99% completion, <2% redundancy. (out of 579 expected single-copy marker genes, 563 were present once, and 11 were present twice). I am especially concerned about the missing genes.
I can think of several possibilities:
a) Limits of CheckM. This software uses a sample of 82 genomes to select marker genes for this phylum: however, there are no available genomes in the same genus as the one I'm working with. It might be possible that this is natural variation.
b) The "single-copy" genes present twice might be due to misassembly. But while the sample was not axenic, the organism of interest was majoritary, and long reads were used as well, so this shouldn't be an issue?
c) The genes were somehow missed by sequencing. (but the depth is pretty good?)
I don't know how likely a) is, so it's difficult for me to assess whether what I have is a normal result, or something that requires additional sequencing in order to make sense. I would appreciate your thoughts!

See this table from the CheckM table to put this number in context. In general, your numbers sound excellent, even compared to isolate genomes.

I'd recommend for context taking some finished microbial genomes (E. coli K12 etc) and running CheckM. I bet you'd be surprised what you find.
I think that you can probably not worry about (c) as a possibility. If you are very worried about (b), you could map reads against your assembly, filter out the ones that don't map well, and assemble those to see what you find. Long reads suggest that your genome is in good shape. You can of course do manual curation of your genome using tools like Anvi'o, or by scaffolding/gap-filling your assembly with long reads or an orthogonal strategy like optical mapping or Hi-C.
I think that the influence of (a) is significant. Overall it sounds like you have a reasonable genome.

Don't read into the CheckM result too much. Genes having one copy in some species may have multiple copies in other species and genes present in some species may be absent in others. This may even happen between strains, let alone between species or between genera. It is rare to see 100% checkM completeness and 0% contamination. A common standard for "near complete" assembly is >90% completeness and <5% contamination. Your assembly is far above this standard.
Meanwhile, beware that checkM only evaluates single-copy genes but multi-copy genes are usually harder to assemble. There might be complex multi-copy genes in the three assembly gaps. Without getting a circular contig, it is hard to say for sure your assembly is really complete.


text:
Compare mapped reads from BWA -MEM and STAR from .bam files
I want to find and compare the results from STAR and BWA- MEM mapped. I have 150bp paired end reads in sorted.bam files in each case and i want to find in each case uniquely mapped reads, number of reads mapped to multiple loci and number of unmappped reads.
Especially for uniquely mapped reads i tried different solutions from biostars and this Obtaining uniquely mapped reads from BWA mem alignment but i dont know what to trust,
Could you please suggest the best way to find mapping results from bam paired end files ?

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.

It would be helpful to describe the genome and a bit of project background regarding what you have already tried. Biostars is a different site to Stackexchange.

I have paired end mapped reads from human in sorted bam files and i want to find the unique mapped , the multimapped and the unmapped reads from BWA and from STAR. I already have the .bam files from each mapper. i Tried samtools flagstat and the solution above but i am not sure if this is correct especially for unique reads

Bam files contain a MAPQ field that informs on whether an alignment is unique or not. With STAR, a MAPQ value of 255 means the read is mapped uniquely and a score lower than 255 means not unique.
You can use samtools on a bam file to count or output the uniquely or multimapped reads by filtering for the MAPQ as described here: https://groups.google.com/forum/#!topic/rna-star/ipHc9FrgaPU
For example, count uniquely mapped reads:
samtools view -c -q 255 myBamFile.bam
Output uniquely mapped reads:
samtools view -q 255 myBamFile.bam
There are differences between STAR and BWA MAPQ output though. see here for details: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/


text:
Chosing an imputation panel for SNP-Chip data?
I have about 1,000 SNP-Chip data (samples) that I'd like to impute over (for the purpose of having more rsIDs to match against GWAS data).
However, I don't know the ancestry of each sample / the ancestry hasn't been recorded in a reliable way.
Is there a 'quick and dirty' method to decide which imputation panel to use based on the genotype data itself? e.g.

sample 1 'looks' British in England and Scotland (GBR),
sample 2 'looks' Colombian in Medellin, Colombia (CLM),
sample 3 'looks' Bengali in Bangladesh (BEB), etc.

Once I know which group each sample matches most closely I can then perform imputation with the appropriate imputation panel.

If you're using an hmm-based method such as for example impute2 or later, then the imputation method will perform the haplotype matching for you. No need to cut down the reference panel. Try to get a big reference panel such as 1000 genomes phase 3 or the haplotype reference consortium (disclaimer: I am an author on both those papers).
But do be careful about relative power: a sample with haplotypes that are better represented in the ref panel than another sample will have better-imputed genotypes than the other sample.

Nice, I didn't know that impute2 did that. Very cool. I take your point that some samples will have better (more extensive and reliable) imputation than others.

If you do want to use the haplotype reference consortium (HRC) panel, then you may need to upload your SNP-chip data to an imputation server, as the HRC panel is partially restricted due to the wishes of the participants.

Also, there has been some recent work on improving imputation speed using very large reference panels: https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1009049

can you recommend an imputation server to give me access to HRC and uses HMM-based method such as for example impute2? You seem to be implying that HRC is the largest reference panel out there, but perhaps that's not true... (sorry for my basic questions).

I cannot vouch for these servers, but they are run by reputable groups (my coauthors on the HRC paper in fact): https://www.sanger.ac.uk/tool/sanger-imputation-service/ and https://imputationserver.sph.umich.edu/index.html#!

If I were you, I would use ADMIXTURE in supervised mode --supervised and use some of the 1000 genomes populations as reference populations. This is a fast and accurate way to obtain broad-scale ancestry proportion estimates for each of your individuals.

Thanks, that would answer my question, but it's a bit of an XY question... I want to get ancestry to do imputation... It looks like there are better ways to do imputation than first selecting ancestry :-)


text:
Can Primer-BLAST be run locally?
Primer-BLAST is a tool for designing target-specific primers, available through the Primer-BLAST web interface.
Is it possible to run the tool locally instead of accessing it via the internet?

Primer-BLAST obviously uses NCBI Blast which we all know about, but it is a pipeline running through Primer3. This is an extremely popular and reliable primer design software, but only does single stranded DNA/RNA. There is a prettier version, Primer3plus here.

To strictly answer the question its not possible because it is a pipeline based on Primer3 which is not written by NCBI, therefore they would not be in a position to distribute it. Nor would Primer3 be able to distribute Blast because that is NCBI. It is simply a web-server.
NCBI and Primer3 would need to come to a joint agreement, which they must have done to present the webserver (at least it would be polite). You can of course replicate the pipeline yourself.
Brief details are here,
https://www.ncbi.nlm.nih.gov/tools/primer-blast/primerinfo.html

It can be done via Docker. This exact pipeline hasn't been distributed via Docker but Blast is available as a Docker image.

New Answer
Just to follow this up, primers4clades is a Docker image aimed at very similar objectives to Primer-Blast. I have not used it, but its worth considering if standalone is essential, e.g. for ease of building a pipeline.
Here's the blurb

Primers4clades provides a fully automated pipeline for the design of PCR primers for cross-species amplification of novel sequences from metagenomic DNA or from uncharacterized organisms belonging to user-specified phylogenetic lineages. It implements an extended CODEHOP strategy based on both DNA and protein multiple alignments of coding genes and evaluates thermodynamic properties of the oligonucleotide pairs, as well as the phylogenetic information content of predicted amplicons, computed from the branch support values of maximum likelihood phylogenies.

I don't know the user, but Docker is a safe running environment that is isolated from the system. There's no harm in checking.
Finally both Blast and Primer3 are available via conda, but not Primer-Blast. In theory, I believe Primer-Blast could be distributed via conda under the rules of conda. In fact, I think, but could be wrong, that anyone could rebuild this pipeline and post it back to conda. @merv will know exactly whether this is 'permitable' and is very likely to correct me if I'm wrong.


text:
Error in installing vcf2maf
I tried to install vcf2maf following this tutorial from GitHub: https://github.com/mskcc/vcf2maf.
I get a repeat error when I try to test out vcf2maf using the code:
perl vcf2maf.pl --input-vcf tests/test.vcf --output-maf tests/test.vep.maf

The error message looks like this:
'''STATUS: Running VEP and writing to: tests/test.vep.vcf
ERROR: Cannot find VEP script under: /Users/user/miniconda3/bin'''

Basically I followed every steps provided on the tutorial.
Does anyone could help me find out what is going wrong?
I am working on macOS Monterey 12.2.1 with an Intel Core i5 processor.

Could you kindly provide details of your hardware and system?

I have the same problem.
I tried conda installation as indicated, however, the error is not solved.
Perl version is v5.36.0.

Install it via anaconda/conda either,
conda install -c bioconda vcf2maf
conda install -c bioconda/label/cf201901 vcf2maf

I would be surprised if that installation fails. A possible issue is Perl version and system version.

Thanks! I think that works!

Cool, please remember to 'accept' and if possible upvote


text:
Affinity maps and forcefields
I know that before performing docking with AutoDock Vina, we have to compute "vina maps" like below
v.set_ligand_from_file('1iep_ligand.pdbqt')
v.compute_vina_maps(center=[15.190, 53.903, 16.917], box_size=[20, 20, 20])

but I do not quite understand what those maps exactly are.
I found from the manual that

...and then compute the affinity maps for each ligand atom types accroding to the Vina forcefield.

but again, what are affinity maps and Vina forcefield. I checked the source code to see that it creates some kind of a grid which will be then used for scoring. What kind of information does this grid data type contains? I would greatly appreciate any help.

Docking is a computational calculation in the field of molecular mechanisms.
It can be argued that whereas it is true one can 'wing it' without knowing the basics, it is best to know them —akin to the fact that it is possible to drive a car fine without knowing the road code or how a car works: many do for sure but it is not a good idea as something will go wrong.
So I would suggest having a gander at the very least at:

https://en.wikipedia.org/wiki/Molecular_mechanics
https://en.wikipedia.org/wiki/Force_field_(chemistry)
https://en.wikipedia.org/wiki/Lennard-Jones_potential
https://en.wikipedia.org/wiki/Solvent_model
https://manual.gromacs.org/documentation/current/reference-manual/topologies/particle-type.html#atom-types

So a forcefield is a sum of different terms that mathematically describes the forces acting between atoms. A thing to note is that these terms may be have some basis in chemistry or some may be just fudge factors to make the data fit.
Regarding the atom types, in a full atom model, each atom is represented as a particle (or node in graph theory speak) and the bonds are fixed (edges).
An element behaves different depending on its hybridisation and charge etc. Therefore, each atom is classified not as an element, but as an atom type. Say, nitrogen in a ring and amine are generally different atom types. With parameterising a small molecule an atomtype may not be a perfect fit.


text:
Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
I implemented Hirschberg's algorithm in python and used the wiki example to verify correct implementation given the scoring parameters and sequences:
Let X = AGTACGCA
    Y = TATGC
Del(X) = -2
Ins(Y) = -2
Sub(x,y) = -1 for mismatch +2 for match

This aligns to:

AGTACGCA
--TATGC-

Which matches the wiki output. I also printed the recursion and it matched the wiki recursive partitioning too. Then I tried it on a harder problem and got a terrible result. It was a mitochondrial Neanderthal DNA sample of about 500 nucleotides and a complete human DNA mitochondrial sample of 16504 nucleotides. The output of trying it on BlastN is here:
https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_Query_40035
My results look something like this where the top is neanderthal and the bottom is the full human genome. It appears though it just matches each character as soon as it finds on.
---C-CA-----A--------------------G-----T----AT---T--G--A----C-TC--------A--C-C-C-AT--CA-----A--C-----A-----A-C-C-----G-C-CA-T----GT-----ATT--T-C------------------G----TAC----A-T-T-AC---------T----------G------------------C--------CA--------------G-----C--CAC--C-A-T---------GA-AT-AT---------T-----------------------G-T-------ACAG---T-A--C-CAT-------AA--------------------T----------TA-C----TT----------------G-----A--C--T---A-C---C----------T--------------------------------------------------------G----T---A---A---------------T---A-C---ATA-----AA--AACC------------------------T--A---A--T--CC-----A---CA-T-CA---A-----------C---C-C-C--C-CCCC-----C-------------CC-A---T--------GCT-T-A---------CA-A-GCAAGCA-C--------AG------CA------A--TCA--AC---C--------------T--------------T-CA----A---C--T--G--T--C-A---------------------T-A----C----A----T---C-AA-------CTA--C-A-ACT--CC--A---------AA----G----A-C-ACC-C--T-------T-A-C--A---C--------CC-----A---------------C------T----A-----G-------------GA-T---A------TC-A----ACA-AA------C--C-----T------A-C----C----CAC-C-----CT-----------TG--A-------C---A----G-T-A-C---A----T-A---G------CA-CA----T-------A-AA------G---------T-----C--A-------T-----T--T--A---CC-----G-------T-------A--C-ATA------G--CA-C---A-----T-T-----A---TA------G---TC---A--AA----T--------C--C--------C-----T-------T-----C-T-----C--G------CCC---------C----------C-A--T-----------G--------GAT-G-----A-----C--------C--------------C------C-------------C---------C-------T----C----A------G------A---T----A-G-----------GG---------------------G--T--------------C----C----C-T--T--G---A---------------------------------------------------------------------------------------------------------------
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAATGTCTGCACAGCCACTTTCCACACAGACATCATAACAAAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCAATACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCAAGCACGCAGCAATGCAGCTCAAAACGCTTAGCCTAGCCACACCCCCACGGGAAACAGCAGTGATTAACCTTTAGCAATAAACGAAAGTTTAACTAAGCTATACTAACCCCAGGGTTGGTCAATTTCGTGCCAGCCACCGCGGTCACACGATTAACCCAAGTCAATAGAAGCCGGCGTAAAGAGTGTTTTAGATCACCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGGCTTTAACATATCTGAACACACAATAGCTAAGACCCAAACTGGGATTAGATACCCCACTATGCTTAGCCCTAAACCTCAACAGTTAAATCAACAAAACTGCTCGCCAGAACACTACGAGCCACAGCTTAAAACTCAAAGGACCTGGCGGTGCTTCATATCCCTCTAGAGGAGCCTGTTCTGTAATCGATAAACCCCGATCAACCTCACCACCTCTTGCTCAGCCTATATACCGCCATCTTCAGCAAACCCTGATGAAGGCTACAAAGTAAGCGCAAGTACCCACGTAAAGACGTTAGGTCAAGGTGTAGCCCATGAGGTGGCAAGAAATGGGCTACATTTTCTACCCCAGAAAACTACGATAGCCCTTATGAAACTTAAGGGTCGAAGGTGGATTTAGCAGTAAACTAAGAGTAGAGTGCTTAGTTGAACAGGGCCCTGAAGCGCGTACACACCGCCCGTCACCCTCCTCAAGTATACTTCAAAGGACATTTAACTAAAACCCCTACGCATTTATATAGAGGAGACAAGTCGTAACATGGTAAGTGTACTGGAAAGTGCACTTGGACGAACCAGAGTGTAGCTTAACACAAAGCACCCAACTTACACTTAGGAGATTTCAACTTAACTTGACCGCTCTGAGCTAAACCTAGCCCCAAAC

The neanderthal sample matches at the end of the human mitochondrial dna. I suspect the gap score penalty is far too costly to insert enough gaps to reach the end of the alignment but I set blast to a linear gap penalty. How do I overcome this? Here is the implementation code for my program:
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

# NWScore
def NWScore(X, Y):
    # Initialize NM Scoring Matrix. Array needs to be the length of X +1 and Y + 1
    score = [[0 for i in range(len(Y)+1)] for j in range(len(X)+1)]
    for j in range(1,len(Y)+1):
        # Insertion scoring penalty column
        score[0][j] = score[0][j-1] + ins("a")
    for i in range(1,len(X)+1):

        # Delete Penalty Row
        score[i][0] = score[i-1][0] + dele("a")
    for i in range(1,len(X)+1):
        for j in range(1,len(Y)+1):
            # Scoring Matrix
            scoreSub = score[i-1][j - 1] + sub(X[i-1], Y[j-1])
            scoreDel = score[i-1][j] + dele(X[i-1])
            scoreIns = score[i][j - 1] + ins(Y[j-1])
            score[i][j] = max(scoreSub, scoreDel, scoreIns)
    if toPrint:
        print(X,Y)
        print("Score")
        for row in score:
            print(row)
    lastLine = score[-1]
    return lastLine

# Hirschberg
def Hirschberg(X, Y):
    Z = ""
    W = ""
    # If sequence X has no more char in the seq
    if len(X) == 0:
        for i in range(len(Y)):
            # Concatenate a gap(?) to Z and Y[i] to W
            Z = Z + '-'
            W = W + Y[i]
    # If sequence Y has no more char in the seq
    elif len(Y) == 0:
        for i in range(len(X)):
            Z = Z + X[i]
            W = W + '-'
    elif len(X) == 1 or len(Y) == 1:
        #If one seq has a single use full global alignment. For now there is a biopython implementation of global alignment

        a = pairwise2.align.globalms(X,Y,2,-1,-2,-2)
        Z = a[0].seqA
        W = a[0].seqB

    else:
        xlen = int(len(X))
        xmid = round(len(X)/2)
        ylen = int(len(Y))
        # Ok so this is the meat of the algorithm. The last few conditionals dealt with the base cases.
        # Scoreleft takes NWScore for the first half of sequence X, and Sequence Y

        ScoreL = NWScore(X[:xmid], Y)
        # I imagine this is to start from the right side.
        ScoreR = NWScore(X[xmid:][::-1], Y[::-1])[::-1]

        ymid = numpy.argmax([x+y for (x,y) in zip(ScoreL,ScoreR)])

        # Recursion. This is divide and conquer. Creating partitions that somehow dynamically reassemble even
        # though portions of it are getting reversed?!
        Z = Z+ Hirschberg(X[:xmid],Y[:ymid])[0]+Hirschberg(X[xmid:xlen],Y[ymid:ylen])[0]
        W = W+Hirschberg(X[:xmid],Y[:ymid])[1] +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])[1]

    return Z,W
seq1 = SeqIO.read("seq1.fasta","fasta")
seq2 = SeqIO.read("seq2.fasta","fasta")

Z,W = Hirschberg(seq1.seq,seq2.seq)
print(Z)
print(W)
score = 0
for i in range(len(Z)):
    if Z[i] == '-' or W[i] == '-':
        score+= insertionPenalty
    elif Z[i]  == W[i]:
        score+=match
    else:
        score+=mismatch

print("Score: ",score)

The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?

I believe so. The fasta file states that it is a complete human mitochondrial dna sequence.

Okay thanks, aligning via amino acids will not work then (which is what I woudl have done). The coding looks fine and its a nicely focused question.

Comparisons First of all I would simply say pipe it through Muscle and then compare the results.
I forget the algorithm advances in the 'new' (now quite old) generation of aligners, but it is worth reading up on the latest Clustal omega, possibly MAFFT and how they view Muscle. I do remember is 2018 Clustal Omega upgrade eulogising the historic innovation of Muscle around accuracy and describing how this was implemented. If that works for you, you simply replicate it in your code.
The modern development in aligning algorithms is how they compare in speed versus number of sequences, but accuracy is always measured.
Amino acids -> nucleotides The reason I don't follow the 'accuracy' issues is because its easier to align via amino acids and then switch back to codons. If you are sticking with this gap, penalty system using amino acids is a good approach and then simply back-translating, but if its a UTR then that will not work.
Brute-forcing the gap penalties is certainly an option (which is the advantages of coding it). Proving why it has worked might be a bit tricky, but its doable. You could simply fire parameters onto the data and use regex (import re) or straight character matching == to assess the improvements. Regex is better to assess improvements, whilst == approach is much faster. It depends on the parameter space (I think regex). Moving to an coded 'optical analysis' is an easy modification and could be useful.

The implication from the comment is there's suboptimal code ...

The raw default output to Biopython's pairwise should be used to debug this. Then parameters changed stepwise to understand the anomalous output.

Skimming the code I don't spot anything obvious. What I don't understand is why the code has been duplicated
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

This is present twice. I think Python would complain about the duplication of functions.
The duplication of variables is fine and probably imports.

A looked through the code a final time and yes there are bugs
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

That can't be right because whatever you call the functions (which you do), each time the output is -3
No matter what you input the output is always -3. If this was OOP and these were class methods in some sort of method chain that  might be different. These are functions which are operating not different to calling the variable direct. Thus if you substitute
... + dele("a")

with
... + deletePenalty

or even
.... + in("a") # I know this isn't supposed to be here but its to make a point

The answer in all cases is -3.
If you asking for bugs, thats got to be involved.
Out of interest are you using an IDE with breakpoints etc...?

I know there are better ways to do it but I successfully used blastn to align the sequences which operates using scoring and global alignment. So it should be possible. I believe I need to create a system to allow for gaps at the beginning and the end of the smaller sequence.

Okay, the implication is there is a bug in the code, I had a brief look.

My mistake, I removed the duplicated code.

Okay, I looked through the code again and in and dele would a key source of bugs

Large sequences of dissimilar size? Sounds like you want a local alignment algorithm. Hirschberg is a global alignment algorithm, which will often perform unexpectedly on large sequences of dissimilar size.


text:
Featureplot by co-expression of some genes
I want to identify MDSCs cells which don't not express HLA-DRs (HLA-DR negative) but express CD14, CD33, and ITGAM
For this, I firstly plotted a feature plot to locate HLA-DR by
cd.features <- list(grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE))

pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')

FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

Which returned

For CD14, CD33, and ITGAM I done
> cd.features=list(c("CD14", "CD33",  "ITGAM"))
> pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')
> FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

>

Which returned

By looking at these two plots I really can not say which part of plot is HLA-DR negative but expressing CD14, CD33, and ITGAM
Do you know a way to do this?
Plotting options don't help because looking at the below plots you can not say which cells are HLA negative and CD14 positive

Will the user who flagged for this post to be closed please add a comment stating why they think this post should be closed? I don't see any problem with it.

My small comment, I'd simply use ggplot2 if I was using R (not hugely helpful). It is a lot more powerful though. @RamRS point: It is a very clear question, nicely presented on a good subject. I can make a good educated guess why someone did this, but its really for meta. Closing questions is a current topic there. It is unlikely the question will be closed - its just a vote. (For the record I upvoted.)

@haci will know this, but will have to get into cellular immunology.

Thinking about it you really want a third axis and a good viewing angle. Plotly might be a better option than featurePlot. Interesting but not my thing.

Very old question, but recently modified, which means that someone might still be looking into this.
I usually prefer to use Ucell as the signatures are calculated within cells, which makes comparisons a bit easier. In addition is allows the formulation of signatures with both positive and negative genes.
hladr_genes <- grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE)
signatures_list <- list(MDSCs = c("CD14+","CD33+","ITGAM+",paste0(hladr_genes,"-")))

> signatures_list
$MDSCs
[1] "CD14+"     "CD33+"     "ITGAM+"    "HLA-DRA-"  "HLA-DRB5-" "HLA-DRB1-"

# calculate signature
pbmc_small <- AddModuleScore_UCell(pbmc_small,signatures_list,ncores = 8,name = "")
# perform smoothing
pbmc_small <- SmoothKNN(pbmc_small, signature.names = names(signatures_list), reduction="pca")

You need to make a transformed variable to the gene sets of interest. I would rescale the data on a range of 0-1. Then you can make a log2 ratio of the sum(set1)/sum(set2). You could explore other transformations, but this will give contrasting values for high set1 and high set2, but neutral if the expression is similar between the sets.

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.


text:
LD Score Regression Derivation hard to follow
I am trying to understand the derivations from Sullivan et al. (2015) in the Supplementary Material. There, it is mentioned in the first page that the least squares estimate of the j-th SNP effect, considering the polygenicity linear model $φ=Xβ + ε$, is $\hat{β_j} := X^T_j φ/N$ . Normally, the least squares solution of $\hat{\beta}$ is $(X^TX)^{-1}X^T\phi$. By comparing the two equations with each other, this means that, in this case, $X^TX=N I_{M\times M}$. However, LD is assumed, therefore $E(X_{ij}X_{ik})\neq 0$, so the covariance matrix cannot be diagonal. Could someone explain to me how the derivation for the beta estimate was performed, please? What am I missing? I checked that the same thing is mentioned in more than one sources such as here.

I still don't understand what you're asking. If you want to try and rephrase your question I can give it a shot, but I'm not sure where the confusion lies right now.

I edited the question, I hope now it is clearer what is the mathematical issue I am having trouble understanding.

Not entirely clear on what you're asking but the beta estimates in LDSC come from GWAS summary statistics. As I understand it, they aren't "derived" so much as just set equivalent to what they linear model defines them as, which is their additive effects on the phenotype as implied by the model. So for each SNP j, it's effect size (beta) is whatever difference one copy of the allele makes on the phenotype (phi), divided by N, the number of SNPs.

Sorry for being late, but what do you mean by the phrasing "equivalent to what the linear model defines them as"? If only the linear model solution was considered, then the estimated beta-s would be the pseudo-inverse of X times X' times phi, namely the least squares solution.

Accepting your answer, I think I misunderstood the phrasing in the supplementary note, it is indeed only a way of defining a weighted effect on the explained output, nothing to do with least squares solution...


text:
How to convert a BAM file into a VCF file
I know that there is a lot of methods to do so, but most of them need to have a reference file. (For example, the GATK force me to provide a reference.)
However, I think the BAM file I am working with already embedded the reference information in it. I tried to use the Samtool (mpileup) to do that, but it seems that the vcf file I get has different formats compared with the VCF documentations.
Could anyone give me some advice?
Also, could anyone provide a method to validate the vcf file without using a reference file?
Thanks a lot!
The following is what I got when running:
samtools view -H

@HD VN:1.5  GO:none SO:coordinate
@SQ SN:chr1 LN:248956422
@SQ SN:chr2 LN:242193529
@SQ SN:chr3 LN:198295559
@SQ SN:chr4 LN:190214555
@SQ SN:chr5 LN:181538259
@SQ SN:chr6 LN:170805979
@SQ SN:chr7 LN:159345973
@SQ SN:chr8 LN:145138636
@SQ SN:chr9 LN:138394717
@SQ SN:chr10    LN:133797422
@SQ SN:chr11    LN:135086622
@SQ SN:chr12    LN:133275309
@SQ SN:chr13    LN:114364328
@SQ SN:chr14    LN:107043718
@SQ SN:chr15    LN:101991189
@SQ SN:chr16    LN:90338345
@SQ SN:chr17    LN:83257441
@SQ SN:chr18    LN:80373285
@SQ SN:chr19    LN:58617616
@SQ SN:chr20    LN:64444167
@SQ SN:chr21    LN:46709983
@SQ SN:chr22    LN:50818468
@SQ SN:chrX LN:156040895
@SQ SN:chrY LN:57227415
@SQ SN:chrM LN:16569
@SQ SN:chr1_KI270706v1_random   LN:175055
@SQ SN:chr1_KI270707v1_random   LN:32032
@SQ SN:chr1_KI270708v1_random   LN:127682
@SQ SN:chr1_KI270709v1_random   LN:66860
@SQ SN:chr1_KI270710v1_random   LN:40176
@SQ SN:chr1_KI270711v1_random   LN:42210
@SQ SN:chr1_KI270712v1_random   LN:176043
@SQ SN:chr1_KI270713v1_random   LN:40745
@SQ SN:chr1_KI270714v1_random   LN:41717
@SQ SN:chr2_KI270715v1_random   LN:161471
@SQ SN:chr2_KI270716v1_random   LN:153799
@SQ SN:chr3_GL000221v1_random   LN:155397
@SQ SN:chr4_GL000008v2_random   LN:209709
@SQ SN:chr5_GL000208v1_random   LN:92689
@SQ SN:chr9_KI270717v1_random   LN:40062
@SQ SN:chr9_KI270718v1_random   LN:38054
@SQ SN:chr9_KI270719v1_random   LN:176845
@SQ SN:chr9_KI270720v1_random   LN:39050
@SQ SN:chr11_KI270721v1_random  LN:100316
@SQ SN:chr14_GL000009v2_random  LN:201709
@SQ SN:chr14_GL000225v1_random  LN:211173
@SQ SN:chr14_KI270722v1_random  LN:194050
@SQ SN:chr14_GL000194v1_random  LN:191469
@SQ SN:chr14_KI270723v1_random  LN:38115
@SQ SN:chr14_KI270724v1_random  LN:39555
@SQ SN:chr14_KI270725v1_random  LN:172810
@SQ SN:chr14_KI270726v1_random  LN:43739
@SQ SN:chr15_KI270727v1_random  LN:448248
@SQ SN:chr16_KI270728v1_random  LN:1872759
@SQ SN:chr17_GL000205v2_random  LN:185591
@SQ SN:chr17_KI270729v1_random  LN:280839
@SQ SN:chr17_KI270730v1_random  LN:112551
@SQ SN:chr22_KI270731v1_random  LN:150754
@SQ SN:chr22_KI270732v1_random  LN:41543
@SQ SN:chr22_KI270733v1_random  LN:179772
@SQ SN:chr22_KI270734v1_random  LN:165050
@SQ SN:chr22_KI270735v1_random  LN:42811
@SQ SN:chr22_KI270736v1_random  LN:181920
@SQ SN:chr22_KI270737v1_random  LN:103838
@SQ SN:chr22_KI270738v1_random  LN:99375
@SQ SN:chr22_KI270739v1_random  LN:73985
@SQ SN:chrY_KI270740v1_random   LN:37240
@SQ SN:chrUn_KI270302v1 LN:2274
@SQ SN:chrUn_KI270304v1 LN:2165
@SQ SN:chrUn_KI270303v1 LN:1942
@SQ SN:chrUn_KI270305v1 LN:1472
@SQ SN:chrUn_KI270322v1 LN:21476
@SQ SN:chrUn_KI270320v1 LN:4416
@SQ SN:chrUn_KI270310v1 LN:1201
@SQ SN:chrUn_KI270316v1 LN:1444
@SQ SN:chrUn_KI270315v1 LN:2276
@SQ SN:chrUn_KI270312v1 LN:998
@SQ SN:chrUn_KI270311v1 LN:12399
@SQ SN:chrUn_KI270317v1 LN:37690
@SQ SN:chrUn_KI270412v1 LN:1179
@SQ SN:chrUn_KI270411v1 LN:2646
@SQ SN:chrUn_KI270414v1 LN:2489
@SQ SN:chrUn_KI270419v1 LN:1029
@SQ SN:chrUn_KI270418v1 LN:2145
@SQ SN:chrUn_KI270420v1 LN:2321
@SQ SN:chrUn_KI270424v1 LN:2140
@SQ SN:chrUn_KI270417v1 LN:2043
@SQ SN:chrUn_KI270422v1 LN:1445
@SQ SN:chrUn_KI270423v1 LN:981
@SQ SN:chrUn_KI270425v1 LN:1884
@SQ SN:chrUn_KI270429v1 LN:1361
@SQ SN:chrUn_KI270442v1 LN:392061
@SQ SN:chrUn_KI270466v1 LN:1233
@SQ SN:chrUn_KI270465v1 LN:1774
@SQ SN:chrUn_KI270467v1 LN:3920
@SQ SN:chrUn_KI270435v1 LN:92983
@SQ SN:chrUn_KI270438v1 LN:112505
@SQ SN:chrUn_KI270468v1 LN:4055
@SQ SN:chrUn_KI270510v1 LN:2415
@SQ SN:chrUn_KI270509v1 LN:2318
@SQ SN:chrUn_KI270518v1 LN:2186
@SQ SN:chrUn_KI270508v1 LN:1951
@SQ SN:chrUn_KI270516v1 LN:1300
@SQ SN:chrUn_KI270512v1 LN:22689
@SQ SN:chrUn_KI270519v1 LN:138126
@SQ SN:chrUn_KI270522v1 LN:5674
@SQ SN:chrUn_KI270511v1 LN:8127
@SQ SN:chrUn_KI270515v1 LN:6361
@SQ SN:chrUn_KI270507v1 LN:5353
@SQ SN:chrUn_KI270517v1 LN:3253
@SQ SN:chrUn_KI270529v1 LN:1899
@SQ SN:chrUn_KI270528v1 LN:2983
@SQ SN:chrUn_KI270530v1 LN:2168
@SQ SN:chrUn_KI270539v1 LN:993
@SQ SN:chrUn_KI270538v1 LN:91309
@SQ SN:chrUn_KI270544v1 LN:1202
@SQ SN:chrUn_KI270548v1 LN:1599
@SQ SN:chrUn_KI270583v1 LN:1400
@SQ SN:chrUn_KI270587v1 LN:2969
@SQ SN:chrUn_KI270580v1 LN:1553
@SQ SN:chrUn_KI270581v1 LN:7046
@SQ SN:chrUn_KI270579v1 LN:31033
@SQ SN:chrUn_KI270589v1 LN:44474
@SQ SN:chrUn_KI270590v1 LN:4685
@SQ SN:chrUn_KI270584v1 LN:4513
@SQ SN:chrUn_KI270582v1 LN:6504
@SQ SN:chrUn_KI270588v1 LN:6158
@SQ SN:chrUn_KI270593v1 LN:3041
@SQ SN:chrUn_KI270591v1 LN:5796
@SQ SN:chrUn_KI270330v1 LN:1652
@SQ SN:chrUn_KI270329v1 LN:1040
@SQ SN:chrUn_KI270334v1 LN:1368
@SQ SN:chrUn_KI270333v1 LN:2699
@SQ SN:chrUn_KI270335v1 LN:1048
@SQ SN:chrUn_KI270338v1 LN:1428
@SQ SN:chrUn_KI270340v1 LN:1428
@SQ SN:chrUn_KI270336v1 LN:1026
@SQ SN:chrUn_KI270337v1 LN:1121
@SQ SN:chrUn_KI270363v1 LN:1803
@SQ SN:chrUn_KI270364v1 LN:2855
@SQ SN:chrUn_KI270362v1 LN:3530
@SQ SN:chrUn_KI270366v1 LN:8320
@SQ SN:chrUn_KI270378v1 LN:1048
@SQ SN:chrUn_KI270379v1 LN:1045
@SQ SN:chrUn_KI270389v1 LN:1298
@SQ SN:chrUn_KI270390v1 LN:2387
@SQ SN:chrUn_KI270387v1 LN:1537
@SQ SN:chrUn_KI270395v1 LN:1143
@SQ SN:chrUn_KI270396v1 LN:1880
@SQ SN:chrUn_KI270388v1 LN:1216
@SQ SN:chrUn_KI270394v1 LN:970
@SQ SN:chrUn_KI270386v1 LN:1788
@SQ SN:chrUn_KI270391v1 LN:1484
@SQ SN:chrUn_KI270383v1 LN:1750
@SQ SN:chrUn_KI270393v1 LN:1308
@SQ SN:chrUn_KI270384v1 LN:1658
@SQ SN:chrUn_KI270392v1 LN:971
@SQ SN:chrUn_KI270381v1 LN:1930
@SQ SN:chrUn_KI270385v1 LN:990
@SQ SN:chrUn_KI270382v1 LN:4215
@SQ SN:chrUn_KI270376v1 LN:1136
@SQ SN:chrUn_KI270374v1 LN:2656
@SQ SN:chrUn_KI270372v1 LN:1650
@SQ SN:chrUn_KI270373v1 LN:1451
@SQ SN:chrUn_KI270375v1 LN:2378
@SQ SN:chrUn_KI270371v1 LN:2805
@SQ SN:chrUn_KI270448v1 LN:7992
@SQ SN:chrUn_KI270521v1 LN:7642
@SQ SN:chrUn_GL000195v1 LN:182896
@SQ SN:chrUn_GL000219v1 LN:179198
@SQ SN:chrUn_GL000220v1 LN:161802
@SQ SN:chrUn_GL000224v1 LN:179693
@SQ SN:chrUn_KI270741v1 LN:157432
@SQ SN:chrUn_GL000226v1 LN:15008
@SQ SN:chrUn_GL000213v1 LN:164239
@SQ SN:chrUn_KI270743v1 LN:210658
@SQ SN:chrUn_KI270744v1 LN:168472
@SQ SN:chrUn_KI270745v1 LN:41891
@SQ SN:chrUn_KI270746v1 LN:66486
@SQ SN:chrUn_KI270747v1 LN:198735
@SQ SN:chrUn_KI270748v1 LN:93321
@SQ SN:chrUn_KI270749v1 LN:158759
@SQ SN:chrUn_KI270750v1 LN:148850
@SQ SN:chrUn_KI270751v1 LN:150742
@SQ SN:chrUn_KI270752v1 LN:27745
@SQ SN:chrUn_KI270753v1 LN:62944
@SQ SN:chrUn_KI270754v1 LN:40191
@SQ SN:chrUn_KI270755v1 LN:36723
@SQ SN:chrUn_KI270756v1 LN:79590
@SQ SN:chrUn_KI270757v1 LN:71251
@SQ SN:chrUn_GL000214v1 LN:137718
@SQ SN:chrUn_KI270742v1 LN:186739
@SQ SN:chrUn_GL000216v2 LN:176608
@SQ SN:chrUn_GL000218v1 LN:161147
@SQ SN:chrEBV   LN:171823
@RG ID:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1  CN:Novogene LB:Wit_Exon9_DEL    PL:illuminPU:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1    SM:Wit_Exon9_DEL
@PG ID:bwa  PN:bwa  CL:bwa mem -t 8 -M -R @RG\tID:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1\tSM:Wit_Exon9_DEL\tLB:Wit_Exon9_DEL\tPU:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1\tPL:illumina\tCN:Novogene /PUBLIC/source/HW/Disease/database/genome/human/hg38/bwa_index_whole/hg38.fa /TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/QC/Wit_Exon9_DEL/Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1_1.clean.fq.gz /TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/QC/Wit_Exon9_DEL/Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1_2.clean.fq.gz  VN:0.7.17-r1188
@PG ID:sambamba CL:view -S -f bam /dev/stdin    PP:bwa  VN:1.0
@PG ID:MarkDuplicates   VN:2.18.9-SNAPSHOT  CL:MarkDuplicates INPUT=[/TJPROJ6/GB_HUMAN/PJ_AI/service/cromwell/cromwell-executions/ToSomatic/fddd223f-25fd-4533-ae86-14e106c58ec5/call-Mapping/shard-0/Mapping/14f30bea-add3-4ead-b7e1-fcb202ae18f5/call-picardMarkdup/inputs/-714370009/Wit_Exon9_DEL.merge.bam] OUTPUT=/TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/Mapping/Wit_Exon9_DEL/Wit_Exon9_DEL.final.bam METRICS_FILE=metric.txt REMOVE_DUPLICATES=false ASSUME_SORTED=true OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 TMP_DIR=[.] VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=4000000 CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag CLEAR_DT=true ADD_PG_TAG_TO_READS=true DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> MAX_OPTICAL_DUPLICATE_SET_SIZE=300000 VERBOSITY=INFO QUIET=false COMPRESSION_LEVEL=5 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json USE_JDK_DEFLATER=false USE_JDK_INFLATER=false   PN:MarkDuplicates
@PG ID:samtools PN:samtools PP:sambamba VN:1.10 CL:samtools view -H /Users/XXXX/Desktop/data/02.Bam/XXXX.final.bam
@PG ID:samtools.1   PN:samtools PP:MarkDuplicates   VN:1.10 CL:samtools view -H /Users/XXXX/Desktop/data/02.Bam/XXXX.final.bam

The output file I get when using Samtools looks like:

This is generally a bad idea....if you have a reference, use it! otherwise your VCF will be wrong, for example, you'll miss all the HOM_VAR variants.

To go from .bam to .vcf, it is fairly easy using bcftools call:
bcftools mpileup -Ob -o <study.bcf> -f <ref.fa> <sample1.bam>
bcftools call -vmO z -o <study.vcf.gz> <study.bcf>

Taken from htslib website.

Does this method is applicable when I only have BAM files without bcf files?

@scott - the bcf (it’s the same as a vcf) is the output, not the input, so you’re fine to use it.

Hi~~ When I actually tried this method, I get some output like [E::faidx_adjust_position] The sequence "chr1_KI270706v1_random" was not found. Do you know how to fix it? Thanks a lot~~~

Recently, I also had a similar question and got my answer from Biostars. See my post and the answer I received: Germline variant calling from BAM without reference.
Basically, I was also wondering if it is possible to call SNVs/indels from aligned BAM files without reference because reference sequence was already used to produce such files. To this end, it appears some BAM files contain useful (but optional) read-level tags such as MD which, if present, would allow one to even reconstruct the entire reference sequence. For more details, please visit the link above.


text:
Suggestion for python graphics library to plot aligned protein sequences
I'm working towards developing some visualization method for multiple protein sequences and subsequences using PYTHON. I want to  show the alignment of many fragment sequences to the original sequences as well as the position of fragments and the amino acids. The amino acids are to be colored by some physico-chemical properties, something similar but not exactly  as this:

Can somebody suggest some Python libraries suitable for this type of work? The output can be graphics file like png or interactive graphics or html, doesn't matter really for me
Thanks

You might have a look at the biotite package, and more specifically to its biotite.sequence.graphics module: https://www.biotite-python.org/apidoc/biotite.sequence.graphics.html#module-biotite.sequence.graphics.
Disclaimer: I haven't used these functionalities myself, apart from using one of their colour schemes for DNA.


text:
Ligand and protein preparation in Vina docking
I am an undergraduate novice, from CS major, intending to run a virtual-screening run using Autodock Vina. A few questions in this regard:

How to choose the best charge model for the ligand?

By default (p. 2), Vina installs Gasteiger-Marsili (GM) partial charges prior to docking. However, GM, being a conformationally-independent approach, is a poor choice when compared to other conformationally-dependent approaches.
In that case, can I choose B3LYP/6-311G/AIM which is the best-performing EEM-parametrization-set, yet published?

What are the upsides and downsides?

How to choose the best charge model for the protein?

Once again, Vina installs Gasteiger-Marsili (GM) partial charges prior to docking and there exist better models.
Is there any better open-source alternative than MOPAC-PM7?

Vina removes non-polar hydrogens and lone pairs from the ligand as well as protein.

Is there any upside to this operation other than decreasing computational complexity?
If time etc. is not a factor, are docking outputs going to be affected to any significant degree in absence of this modification?

@M__ Narrowed down my previous question

Ideally StackExchange questions should be related to a single problem. Can you please give more details about your project / problem? If this is for an assignment (or similar canned problem), please add the tag 'homework'. If you're interested in general discussion rather than a solution to a specific problem, I recommend you instead post to https://reddit.com/r/bioinformatics

Your opening statement forces me to ask: have you done a re-docking test with vanilla settings and what was wrong with it? Or are you well past that step?
Gasteiger partial charges are more than adequate. The drude particle mod in some FF addresses the polarisability problem, but is uncommon. Some people run QM charge calculations in Gaussian, for better partial charges, but most folk don't. The major issues in docking stem from protein —charged active site sidechains, especially involving cysteine and histidine, loops that ought to have been flipped to close up on the complex etc


text:
LRT test in deseq2 on multiple condition
I have some doubts when i'm running LRT
This is my condition data-frame small subset
 dput(head(metadata,18))
structure(list(FAB = structure(c(1L, 1L, 1L, 1L, 6L, 2L, 6L,
4L, 3L, 5L, 7L, 3L, 2L, 7L, 6L, 2L, 6L, 2L), .Label = c("C",
"M0", "M1", "M2", "M3", "M4", "M5"), class = "factor")), row.names = c("H1",
"H2", "H3", "H4", "TCGA-AB-2856", "TCGA-AB-2849", "TCGA-AB-2971",
"TCGA-AB-2930", "TCGA-AB-2891", "TCGA-AB-2872", "TCGA-AB-2851",
"TCGA-AB-3011", "TCGA-AB-2949", "TCGA-AB-2981", "TCGA-AB-2965",
"TCGA-AB-2822", "TCGA-AB-2828", "TCGA-AB-2959"), class = "data.frame")

Now I m running LRT on this.  In pairwise comparison condition I can set the reference level as such
coldata$condition<- relevel(coldata$condition, ref = "H") where  H` is my control set with whom I compare all other factors.
Now if I have to do the same with LRT is it correct to set the reference level?
Second doubt
If i have to run LRT using different reference level how to do that?
dds <- DESeq(dds, reduced=~condition, test="LRT")
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)


text:
Regression plot of a continuous trait - is there a binary equivalent?
Within the GWAS wikipedia page, you can view the following plot:

"Illustration of a simulated genotype by phenotype regression for a
single SNP. Each dot represents an individual. A GWAS of a continuous
trait essentially consists of repeating this analysis at each SNP."

This is a regression plot showing on the y-axis the phenotype value for a continuous trait. Is it possible to create a similar plot with a binary trait - I can't think what this would look like, or if there is an equivalent plot I could create.

I have not really seen this type of plot used when a binary outcome is being analyzed. I think in part this is due to the fact that the plot will be rather boring and difficult to interpret as there are only a limited set of possible values (4 if the response and independent variable are both binary), and without introducing some level of random noise, all of the data observations with the same X-Y value pair would overlap.
The situation is a little better when the independent variable is continuous: then you can show more of a distribution. In the figure below, pulled from Figure 2 of this paper explaining how to visualize the discrimination abilities of a logistic regression model, the authors have plotted the risk score produced by the logistic regression model for each patient (x-axis) with the patient's actual outcome on the y-axis. Note that some random noise has been introduced to each observation's y-value so that they do not cover up each other. This plot shows that the risk score is able to somewhat separate the two outcomes.

If you are dealing with a categorical independent variable like the figure included in your question, I think it is better to visualize the proportion of the outcome for each independent variable category like this:

Here the outcome is whether the borrower defaulted on their loan (1 = yes, 0 = no) and the independent variable is whether they are a student (yes/no). For the y-axis, I have calculated the proportion of loan recipients in each group that defaulted. This could easily be done with GWAS data such as genotype vs whether the patient died of heart disease.
I also suggest you check out this answer on how to visualize a contingency table, as the data you are describing (assuming the independent variable stays categorical) is a contingency table.

Thank you for the thoughtful answer - that was really helpful!


text:
How to sort legends alphabetically in plotly dash with category_orders function
I am using Plotly dash for creating a scatter plot. The legends I am getting in the final figure are randomly placed. I want to sort labels alphabetically from A to Z with category_orders.
fig = px.box(
        df,
        x=selected_x,
        y=selected_y,
        points='all',
        hover_data=hover_data,
        color=colour_by,
        width=800,
        height=600,
        labels={
            selected_y: "{} {}".format(selected_y_gene, selected_y),
            selected_x: "{} {}".format(selected_x_gene, selected_x),
        },

Before you construct the figure could you try sorting the values in the dataframe in the order you require? Thus set the True/False flag according to what you want.
orderedDF = df.sort_values(
  by=col_label,
  ascending=True)

then add the following into the plot denoting the labels:
sort=False

Thus added below.
fig = px.box(
    orderedDF,
    x=selected_x,
    y=selected_y,
    points='all',
    hover_data=hover_data,
    color=colour_by,
    width=800,
    height=600,
    labels={
        selected_y: "{} {}".format(selected_y_gene, selected_y),
        selected_x: "{} {}".format(selected_x_gene, selected_x),
# I think it should go here
        sort=False
    }
    )

The gist is you are wanting to stop the automated re-ordering and force the dataframe order onto it.


text:
What is the best tools to find all the somatic mutations?
I am a beginner in this field. I would like to know which is the best way to get all the mutations from a certain sample. At this point, I am considering using GATK (here) and Maftools (here). I don't know which method is better, or are there any other better tools you would recommend?

MAFtools is not a variant calling tool. Google "benchmark somatic variant callers" and read some papers - MuTect2, SomaticSniper, VarScan2, Strelka etc are some relevant tools.

What organism are you working on, i.e. your sample? This has an enormous impact on the approach and the tools.


text:
Pulling out genes in a scatterplot
I'm comparing gene expression among 2 different datasets (in vivo and in vitro) I have made a heatmap showing the correlation for each and then plotted the data frame to create a scatterplot. Now I want to just pull out KP4 from each data set and do the same scatterplot but then pull put genes that aren't on the linear line. How do I do that?
###how to compare expression across both datasets??
# Plot the correlation heatmap of intersect genes between two datasets
dat_human.df_num <- apply(dat_human.df,2,as.numeric)
rownames(dat_human.df_num) <- rownames(dat_human.df)
colnames(dat_human.df_num) <- c("KP4_invitro", "L3.3_invitro", "T3M4_invitro")
colnames(countmatrix_clean_cl) <- c("KP4_invivo", "T3M4_invivo", "L3.3_invivo")
intersect.genes <- intersect(rownames(countmatrix_clean_cl), rownames(dat_human.df_num))
length(intersect.genes) # 18786
#test <- data.frame(cbind(countmatrix_clean_cl[intersect.genes,], dat_human.df_num[intersect.genes,]))
combined_data <- data.frame(cbind(countmatrix_clean_cl[intersect.genes,], dat_human.df_num[intersect.genes,]))
head(combined_data)

melted_corr_mat <- melt(cor(combined_data))
head(melted_corr_mat)

ggplot(data=melted_corr_mat, aes(x=Var1, y=Var2, fill=value)) +
  geom_tile(color="white") +
  scale_fill_gradient2(low="blue", high="red", mid="white", midpoint = mean(melted_corr_mat$value), space="Lab") +
  geom_text(aes(Var2, Var1, label = round(value, digits=2)),
            color = "black", size = 4) +
  theme_minimal() +
  labs(x="",y="") +
  theme(axis.text.x=element_text(angle=45, vjust=1, hjust=1)) +
  coord_fixed()

### scatterplot for each cell line in each condition (in-vivo vs. in-vitro)
head(combined_data)
         KP4_invivo T3M4_invivo L3.3_invivo KP4_invitro L3.3_invitro T3M4_invitro
TSPAN6       24.084      12.762      24.877      31.043       69.536       30.780
TNMD          0.000       0.000       0.000       0.000        0.000        0.000
DPM1         56.015     100.402      60.142     101.408      149.830      214.116
SCYL3         3.470       5.474       8.940       6.889        9.599       10.494
C1orf112      7.910      15.051      12.932      19.363       28.695       38.230
FGR           0.042       0.776       0.128       0.010        0.392        0.361

plot(combined_data)

Why don't you just perform differential expression, that tells you which genes are different, hence not on a diagonal?

I'm not sure how to perform differential expression on bulk?


text:
Help using BioPython's IC_Chain functionality for getting protein 3D coordinates from torsion angles
I have been successfully using BioPython's functionality for extracting protein torsion angles using the IC_Chain module. However, I am currently looking to use the inverse functionality as follows:
from Bio.PDB import PDBParser, internal_coords

# Read PDB entry.
parser = PDBParser(PERMISSIVE=1, QUIET=1)
structure = parser.get_structure('', 'example.pdb')

# Generate protein internal coordinates (angles).
ic_chain_bound = internal_coords.IC_Chain(structure)
ic_chain_bound.atom_to_internal_coordinates()

# Generate back 3D coordinates from angles.
ic_chain_bound.internal_to_atom_coordinates()

but this fails with the following error:
Traceback (most recent call last):
  File "test_biopython_internal_coords.py", line 97, in <module>
    main(args)
  File "test_biopython_internal_coords.py", line 89, in main
    ic_chain_bound.internal_to_atom_coordinates()
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 661, in internal_to_atom_coordinates
    verbose=verbose, start=start, fin=fin
  File "software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 415, in assemble_residues
    Dict[AtomKey, numpy.array], ric.assemble(verbose=verbose)
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 1732, in assemble
    atomCoords = self.get_startpos()
  File "/software/anaconda/installation/envs/alphafold/lib/python3.7/site-packages/Bio/PDB/internal_coords.py", line 1632, in get_startpos
    sp = self.residue.parent.internal_coord.initNCaC.get(self.rbase, None)
AttributeError: 'NoneType' object has no attribute 'initNCaC'

Is this not the intended use? Could someone help me figure out why this is failing?
I'm using this PDB file in my example, but the code previously failed for all other PDB files I tried.
Thanks!

Your ic_chain_bound.residue.parent.internal_coord is None... So your ic_chain_bound.residue.parent needs to be given coordinates... help and dir on it may help

@MatteoFerla set_residues() ---> set_residue() ??? or something else ? trying to figure out from https://biopython.org/docs/latest/api/Bio.PDB.internal_coords.html but I am lost

here https://biopython.org/docs/dev/api/Bio.PDB.internal_coords.html seems better

It turns out that in order to generate the internal coordinates for a protein structure, one shouldn't manually create an IC_Chain object, but rather just use the atom_to_internal_coordinates() method of the structure and the internal_to_atom_coordinates() one for the inverse functionality.
structure.atom_to_internal_coordinates()

# Now you have access to the bond angles and lengths
# of the first chain of structure as follows:
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord

# Dictionary mapping 3 AtomKeys to a Hedron
# (with len12, len23, angle fields):
# ic_chain.hedra

# Dictionary mapping 4 AtomKeys to a Dihedraon
# (with angle field for torsion angles):
# ic_chain.dihedra

structure.internal_to_atom_coordinates()

Maybe you know the answer to this related (and bountied) question?

not able to grasp the question, tried to figure out a way to show the internal_coordinates:
from Bio.PDB import PDBParser, PICIO, PDBIO

# Read PDB entry.
parser = PDBParser(PERMISSIVE=1, QUIET=1)
structure = parser.get_structure('', 'example.pdb')

print(dir(structure)) # 'atom_to_internal_coordinates' and
                      # 'internal_to_atom_coordinates' seems
                      # to be structure methods

# Generate protein internal coordinates (angles).

structure.atom_to_internal_coordinates(verbose = True)

PICIO.write_PIC(structure[0]['A'], 'internal')

# Generate back 3D coordinates from angles.

pippo = PICIO.read_PIC('internal')

pippo.internal_to_atom_coordinates(verbose = True)

print(pippo)

io = PDBIO()

io.set_structure(pippo)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)

compare example.pdb to internal and atom_coord.pdb;
have a look at internal  file contents if you comment out
structure.atom_to_internal_coordinates(verbose = True) (around line 17 of
code, # Generate protein internal coordinates (angles). block-code )