text:
Where to find asymmetric nucleotide substitution matrix with IUPAC encodings?
Update: I submitted a pull request to the Biostrings repo. The functionality I describe in my question and answer can now be implemented with nucleotideSubstitutionMatrix(symmetric = FALSE).

I am using the pairwiseAlignment function from the Biostrings package to calculate the distance between a consensus sequence and a set of other sequences. Some example code:
# R v4.2.0
library(Biostrings)   # v2.63.3
library(BiocGenerics) # v0.41.2
library(tidyverse)    # v1.3.1

# consensus sequence
con <- DNAStringSet(c("GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC"))

# subject sequences
seq <- DNAStringSet(c("ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC",
                      "GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC"))

# substitution matrix
sub <- nucleotideSubstitutionMatrix(match = 0, mismatch = 1)

# get alignments
sapply(X = 1:length(seq),
       FUN = function(i) {
         pairwiseAlignment(pattern = con,
                           subject = seq[i],
                           substitutionMatrix = sub) %>%
         score
      }
    ) %>% as_tibble %>% mutate(sequence = as.character(seq))

The resulting tibble looks like this:
# A tibble: 5 × 2
  value sequence
  <dbl> <chr>
1  -3.5 ATTGTGATTTGCTTTCAAATTAGTATCTTTAAACCATTGAAAACAGC
2  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
3  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
4  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC
5  -0.5 GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC

Something isn't quite right. Take a look at the alignment of con against seq[2].
con:    GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC
        |||||||||||||||| ||||||||||||||||||||||||||||||
seq[2]: GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC

Note that the only difference is at the position in the consensus sequence where there is ambiguous IUPAC encoding. Because R can represent A or G, the score for seq[2] should be 0 (a perfect match). If the ambiguous base was instead in seq[2] and con was unambiguous, then the score penalty would be appropriate. The discrepancy here is due to the symmetry of the substitution matrix sub:
           A          C          G          T          M          R          W          S          Y          K
A  0.0000000 -1.0000000 -1.0000000 -1.0000000 -0.5000000 -0.5000000 -0.5000000 -1.0000000 -1.0000000 -1.0000000
C -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.5000000 -1.0000000
G -1.0000000 -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000
T -1.0000000 -1.0000000 -1.0000000  0.0000000 -1.0000000 -1.0000000 -0.5000000 -1.0000000 -0.5000000 -0.5000000
M -0.5000000 -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -1.0000000
R -0.5000000 -1.0000000 -0.5000000 -1.0000000 -0.7500000 -0.5000000 -0.7500000 -0.7500000 -1.0000000 -0.7500000
W -0.5000000 -1.0000000 -1.0000000 -0.5000000 -0.7500000 -0.7500000 -0.5000000 -1.0000000 -0.7500000 -0.7500000
S -1.0000000 -0.5000000 -0.5000000 -1.0000000 -0.7500000 -0.7500000 -1.0000000 -0.5000000 -0.7500000 -0.7500000
Y -1.0000000 -0.5000000 -1.0000000 -0.5000000 -0.7500000 -1.0000000 -0.7500000 -0.7500000 -0.5000000 -0.7500000
K -1.0000000 -1.0000000 -0.5000000 -0.5000000 -1.0000000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.5000000
V -0.6666667 -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.8333333
H -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.6666667 -0.8333333
D -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.8333333 -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.6666667
B -1.0000000 -0.6666667 -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.8333333 -0.6666667 -0.6666667 -0.6666667
N -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.7500000
           V          H          D          B     N
A -0.6666667 -0.6666667 -0.6666667 -1.0000000 -0.75
C -0.6666667 -0.6666667 -1.0000000 -0.6666667 -0.75
G -0.6666667 -1.0000000 -0.6666667 -0.6666667 -0.75
T -1.0000000 -0.6666667 -0.6666667 -0.6666667 -0.75
M -0.6666667 -0.6666667 -0.8333333 -0.8333333 -0.75
R -0.6666667 -0.8333333 -0.6666667 -0.8333333 -0.75
W -0.8333333 -0.6666667 -0.6666667 -0.8333333 -0.75
S -0.6666667 -0.8333333 -0.8333333 -0.6666667 -0.75
Y -0.8333333 -0.6666667 -0.8333333 -0.6666667 -0.75
K -0.8333333 -0.8333333 -0.6666667 -0.6666667 -0.75
V -0.6666667 -0.7777778 -0.7777778 -0.7777778 -0.75
H -0.7777778 -0.6666667 -0.7777778 -0.7777778 -0.75
D -0.7777778 -0.7777778 -0.6666667 -0.7777778 -0.75
B -0.7777778 -0.7777778 -0.7777778 -0.6666667 -0.75
N -0.7500000 -0.7500000 -0.7500000 -0.7500000 -0.75

As part of a discussion on Github, I provided the code to manually reassign elements of sub to create an asymmetric substitution matrix:
sub["N", "A"] <- 0
sub["N", "T"] <- 0
sub["N", "G"] <- 0
sub["N", "C"] <- 0
sub["M", "A"] <- 0
sub["M", "C"] <- 0
sub["R", "A"] <- 0
sub["R", "G"] <- 0
sub["W", "A"] <- 0
sub["W", "T"] <- 0
sub["S", "C"] <- 0
sub["S", "G"] <- 0
sub["Y", "C"] <- 0
sub["Y", "T"] <- 0
sub["K", "G"] <- 0
sub["K", "T"] <- 0
sub["V", "A"] <- 0
sub["V", "C"] <- 0
sub["V", "G"] <- 0
sub["H", "A"] <- 0
sub["H", "C"] <- 0
sub["H", "T"] <- 0
sub["D", "A"] <- 0
sub["D", "G"] <- 0
sub["D", "T"] <- 0
sub["B", "C"] <- 0
sub["B", "G"] <- 0
sub["B", "T"] <- 0

This covers all cases where the pattern (con) contains ambiguous bases and the subjects (seq) are unambiguous, changing the matrix so the subjects are not penalized. However, these edits do not account for cases where both the pattern and the subject have ambiguous encodings at the same position; e.g. ["M","V"] and ["V","M"] have the same value associated with them, even though V (A || C || G) is a superset of M (A || C). There are also more nuanced cases where the set of bases represented by a pair of ambiguous encodings overlap incompletely in both directions, as with H (A || C || T) and D (A || G || T). Given enough time, I could probably figure out the appropriate scores for each case to complete the asymmetric substitution matrix, but I'm surely reinventing the wheel. Can someone supply, as a standalone table or as part of a codebase, an asymmetric nucleotide substitution matrix that satisfies the criteria outlined herein?

Background
I've dug into the source code for the nucleotideSubstitutionMatrix() function from the Biostrings package. Here is the code for the full function from GitHub:
nucleotideSubstitutionMatrix <- function(match = 1, mismatch = 0, baseOnly = FALSE, type = "DNA")
{
  "%safemult%" <- function(x, y) ifelse(is.infinite(x) & y == 0, 0, x * y)
  type <- match.arg(type, c("DNA", "RNA"))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop("'match' and 'mismatch' must be non-missing numbers")
  if (baseOnly)
    letters <- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters <- IUPAC_CODE_MAP
  if (type == "RNA")
    names(letters) <- chartr("T", "U", names(letters))
  nLetters <- length(letters)
  splitLetters <- strsplit(letters,split="")
  submat <- matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  for(i in 1:nLetters)
    for(j in i:nLetters)
      submat[i,j] <- submat[j,i] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "=="))
  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)
}

Dissecting this, we can see that the matrix is constructed in four main steps.
1. The IUPAC_CODE_MAP named character vector is split into a list of character vectors, where each named element in the list gives an IUPAC character and its constituent DNA bases.
> IUPAC_CODE_MAP
     A      C      G      T      M      R      W      S      Y      K      V      H      D      B      N
   "A"    "C"    "G"    "T"   "AC"   "AG"   "AT"   "CG"   "CT"   "GT"  "ACG"  "ACT"  "AGT"  "CGT" "ACGT"
> splitLetters$A
[1] "A"
> splitLetters$V
[1] "A" "C" "G"
> splitLetters$N
[1] "A" "C" "G" "T"

2. submat is initiated as an $n × n$ matrix of zeros, where $n$ is length(names(IUPAC_CODE_MAP)).
> matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  A C G T M R W S Y K V H D B N
A 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
T 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
M 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
R 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
W 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
S 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Y 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
K 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
V 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
H 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
D 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
B 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
N 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

3. Each element of the matrix is calculated by counting the overlapping elements in the arrays represented by the IUPAC labels given as row and column names.
> for(i in 1:nLetters)
>     for(j in i:nLetters)
>         submat[i,j] <- submat[j,i] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "=="))
> submat
          A         C         G         T         M         R         W         S         Y         K         V         H         D
A 1.0000000 0.0000000 0.0000000 0.0000000 0.5000000 0.5000000 0.5000000 0.0000000 0.0000000 0.0000000 0.3333333 0.3333333 0.3333333
C 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.3333333 0.3333333 0.0000000
G 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.5000000 0.0000000 0.5000000 0.3333333 0.0000000 0.3333333
T 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 0.3333333 0.3333333
M 0.5000000 0.5000000 0.0000000 0.0000000 0.5000000 0.2500000 0.2500000 0.2500000 0.2500000 0.0000000 0.3333333 0.3333333 0.1666667
R 0.5000000 0.0000000 0.5000000 0.0000000 0.2500000 0.5000000 0.2500000 0.2500000 0.0000000 0.2500000 0.3333333 0.1666667 0.3333333
W 0.5000000 0.0000000 0.0000000 0.5000000 0.2500000 0.2500000 0.5000000 0.0000000 0.2500000 0.2500000 0.1666667 0.3333333 0.3333333
S 0.0000000 0.5000000 0.5000000 0.0000000 0.2500000 0.2500000 0.0000000 0.5000000 0.2500000 0.2500000 0.3333333 0.1666667 0.1666667
Y 0.0000000 0.5000000 0.0000000 0.5000000 0.2500000 0.0000000 0.2500000 0.2500000 0.5000000 0.2500000 0.1666667 0.3333333 0.1666667
K 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.2500000 0.2500000 0.2500000 0.2500000 0.5000000 0.1666667 0.1666667 0.3333333
V 0.3333333 0.3333333 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.1666667 0.1666667 0.3333333 0.2222222 0.2222222
H 0.3333333 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.1666667 0.3333333 0.1666667 0.2222222 0.3333333 0.2222222
D 0.3333333 0.0000000 0.3333333 0.3333333 0.1666667 0.3333333 0.3333333 0.1666667 0.1666667 0.3333333 0.2222222 0.2222222 0.3333333
B 0.0000000 0.3333333 0.3333333 0.3333333 0.1666667 0.1666667 0.1666667 0.3333333 0.3333333 0.3333333 0.2222222 0.2222222 0.2222222
N 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000 0.2500000
          B    N
A 0.0000000 0.25
C 0.3333333 0.25
G 0.3333333 0.25
T 0.3333333 0.25
M 0.1666667 0.25
R 0.1666667 0.25
W 0.1666667 0.25
S 0.3333333 0.25
Y 0.3333333 0.25
K 0.3333333 0.25
V 0.2222222 0.25
H 0.2222222 0.25
D 0.2222222 0.25
B 0.3333333 0.25
N 0.2500000 0.25

4. Finally, the base substitution matrix is converted to a scoring matrix given the user-defined values for match and mismatch.

Answer
Concerning step 3, the Biostrings approach asks whether each element of each vector matches each element of the other vector, and then records the proportion of TRUE values in the  matrix of results. Note that the code submat[i,j] <- submat[j,i] creates a symmetric matrix $A$. However, even if $[A]_{ij}$ and $[A]_{ji}$ were generated by separate calls, the commutativity of the mean(outer(i, j, "==")) computation means that $A$ will always be symmetric:
> outer(splitLetters$K, splitLetters$V, "==")
      [,1]  [,2]  [,3]
[1,] FALSE FALSE  TRUE
[2,] FALSE FALSE FALSE
> mean(outer(splitLetters$K, splitLetters$V, "=="))
[1] 0.1666667

> outer(splitLetters$V, splitLetters$K, "==")
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE FALSE
> mean(outer(splitLetters$V, splitLetters$K, "=="))
[1] 0.1666667

Instead of checking for pairwise equality of elements (==), we can check if a given element from the query is contained in the set of elements in the consensus / subject (%in%).
> outer(splitLetters$K, splitLetters$V, "%in%")
      [,1]  [,2]  [,3]
[1,]  TRUE  TRUE  TRUE
[2,] FALSE FALSE FALSE
> mean(outer(splitLetters$K, splitLetters$V, "%in%"))
[1] 0.5

> outer(splitLetters$V, splitLetters$K, "%in%")
      [,1]  [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,]  TRUE  TRUE
> mean(outer(splitLetters$V, splitLetters$K, "%in%"))
[1] 0.3333333

This is not a commutative operation, so the resulting substitution matrix is asymmetric.
asym_NSM <- function(match = 1, mismatch = 0, baseOnly = FALSE, type = "DNA")
{
  "%safemult%" <- function(x, y) ifelse(is.infinite(x) & y == 0, 0, x * y)
  type <- match.arg(type, c("DNA", "RNA"))
  if (!isSingleNumber(match) || !isSingleNumber(mismatch))
    stop("'match' and 'mismatch' must be non-missing numbers")
  if (baseOnly)
    letters <- IUPAC_CODE_MAP[DNA_BASES]
  else
    letters <- IUPAC_CODE_MAP
  if (type == "RNA")
    names(letters) <- chartr("T", "U", names(letters))
  nLetters <- length(letters)
  splitLetters <- strsplit(letters,split="")
  submat <- matrix(0, nrow = nLetters, ncol = nLetters, dimnames = list(names(letters), names(letters)))
  for(i in 1:nLetters) {
    for(j in i:nLetters) {
      submat[i,j] <- mean(outer(splitLetters[[i]], splitLetters[[j]], "%in%"))
      submat[j,i] <- mean(outer(splitLetters[[j]], splitLetters[[i]], "%in%"))
    }
  }

  abs(match) * submat - abs(mismatch) %safemult% (1 - submat)

}

> asym_NSM()
          A         C         G         T         M         R         W         S         Y         K         V         H         D
A 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000
C 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000
G 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000
T 0.0000000 0.0000000 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000 1.0000000 1.0000000 0.0000000 1.0000000 1.0000000
M 0.5000000 0.5000000 0.0000000 0.0000000 1.0000000 0.5000000 0.5000000 0.5000000 0.5000000 0.0000000 1.0000000 1.0000000 0.5000000
R 0.5000000 0.0000000 0.5000000 0.0000000 0.5000000 1.0000000 0.5000000 0.5000000 0.0000000 0.5000000 1.0000000 0.5000000 1.0000000
W 0.5000000 0.0000000 0.0000000 0.5000000 0.5000000 0.5000000 1.0000000 0.0000000 0.5000000 0.5000000 0.5000000 1.0000000 1.0000000
S 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 1.0000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000
Y 0.0000000 0.5000000 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000 1.0000000 0.5000000
K 0.0000000 0.0000000 0.5000000 0.5000000 0.0000000 0.5000000 0.5000000 0.5000000 0.5000000 1.0000000 0.5000000 0.5000000 1.0000000
V 0.3333333 0.3333333 0.3333333 0.0000000 0.6666667 0.6666667 0.3333333 0.6666667 0.3333333 0.3333333 1.0000000 0.6666667 0.6666667
H 0.3333333 0.3333333 0.0000000 0.3333333 0.6666667 0.3333333 0.6666667 0.3333333 0.6666667 0.3333333 0.6666667 1.0000000 0.6666667
D 0.3333333 0.0000000 0.3333333 0.3333333 0.3333333 0.6666667 0.6666667 0.3333333 0.3333333 0.6666667 0.6666667 0.6666667 1.0000000
B 0.0000000 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667 0.6666667
N 0.2500000 0.2500000 0.2500000 0.2500000 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000 0.5000000 0.7500000 0.7500000 0.7500000
          B N
A 0.0000000 1
C 1.0000000 1
G 1.0000000 1
T 1.0000000 1
M 0.5000000 1
R 0.5000000 1
W 0.5000000 1
S 1.0000000 1
Y 1.0000000 1
K 1.0000000 1
V 0.6666667 1
H 0.6666667 1
D 0.6666667 1
B 1.0000000 1
N 0.7500000 1

Additional answer on asymetric IUPAC score matrix.
To get the scores with an asymentric score matrix. This is also implemented in the Bioconductor package MSA2dist (https://www.bioconductor.org/packages/release/bioc/html/MSA2dist.html).
library(Biostrings)
library(MSA2dist)

# consensus sequence and subject sequences
seqs <- setNames(DNAStringSet(c("GTTGTGATTTGCTTTCRAATTAGTATCTTTGAACCATTGAAAACAAC",
"GTTGTGATTTGCTTTCAAATTAGTATCTTTGAACCATTGAAAACAAC")), c("con", "seq"))

# define score matrix
scoreMatrix <- MSA2dist::iupacMatrix()
scoreMatrix["A", "R"]
# [1] 0.5
scoreMatrix["R", "A"] <- 1

MSA2dist::dnastring2dist(seqs, score=scoreMatrix)
MSA2dist::dnastring2dist(rev(seqs), score=scoreMatrix)

# or here using the asymetric score matrix from Biostrings

ansm <- nucleotideSubstitutionMatrix(symmetric = FALSE)

#note: MSA2dist score is normalized by analysed sites
MSA2dist::dnastring2dist(seqs, score=ansm)
MSA2dist::dnastring2dist(rev(seqs), score=ansm)

pairwiseAlignment(seqs[[1]], seqs[[2]], substitutionMatrix = ansm)
pairwiseAlignment(seqs[[2]], seqs[[1]], substitutionMatrix = ansm)

Since MSA2dist::dnastring2dist was until now calculating all possible pairwise combinations from top to bottom (leaving out the already seen combinations), now I will include the option to really calculate all pairwise combinations to take into account asymetric scores and sequence position.


text:
Subsetting out prefix from genes and merging
I'm working with PTX samples (human and mouse), they're aligned to both transcriptomes and so gene names contain hg38- and mm10- at the beginning of each gene. I have eliminated teh prefixes because I want to integrate this dataset with a dataset that does not contain prefixes. I have annotated each cell type accordingly and have subsetted out the mouse cells.
human <- grep(pattern = "^hg38-", x = rownames(x = scData), value = TRUE)
mouse <- grep(pattern = "^mm10-", x = rownames(x = scData), value = TRUE)

mouse_1 <- substring(mouse,6)

human_1 <- substring(human,6)

human_only <- subset(PTX_no_prefix, features = human_1)

mouse_only <- subset(PTX_no_prefix, features = mouse_1)

Idents(mouse_only) <- "cell_ann"
TME_PTX <- subset(mouse_only, idents = c("Fibroblast_mouse", "Myeloid_mouse", "Endothelial_mouse", "Pericytes_mouse", "Neutrophil_mouse"))

I'm now trying to merge two datasets together but receive an error.
combined_het <- merge(TME_het_PTX, y = skin_only_ann)

scData <- combined_het

# remove cells with <200 RNA molecules, or >6000 molecules, or >30% mito
scData <- subset(scData, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & percent.mt_mouse < 30)

#scData <- NormalizeData(scData, normalization.method = "LogNormalize", scale.factor = 10000)
#scData <- FindVariableFeatures(scData, selection.method = "vst", nfeatures = 2000)

all.genes <- rownames(scData)
GLOBALS_MAXSIZE <- 4000
options(future.globals.maxSize = GLOBALS_MAXSIZE * 1024^2)
scData <- SCTransform(scData, vars.to.regress = "nCount_RNA", verbose = FALSE)
#scData <- ScaleData(scData, features = VariableFeatures(object = scData), vars.to.regress = c("nCount_RNA"))
scData <- RunPCA(scData, features = VariableFeatures(object = scData))

combined_het <- merge(TME_PTX, y = skin_only_ann)
Error in intI(i, n = x@Dim[1], dn[[1]], give.dn = FALSE) :
  invalid character indexing

I checked the gene names from each dataset to make sure they're in the same format and they are.
head(rownames(skin_only_ann))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm37686" "Gm1992"  "Gm37329"
head(rownames(TME_PTX))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm37686" "Gm38148" "Gm19938"

possibly related? https://github.com/satijalab/seurat/issues/4452

source: Stack Exchange

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

sort: votes


text:
Bowtie2 mapping with multiple indexes
Background:
We performed NGS using cells collected from mice in a xenotransplantation study.
As such, the FASTQ files contain reads of DNA from both mice and human cells.
I expect ~30% of reads are contaminated by mice cells.
Question:
Is it wrong to directly align using the indexed human reference with bowtie2?
Is there a way to use two indexes at once with bowtie2? The goal here would be to get the % of reads that map to each genome.

I am not aware of a method using two indices in bowtie2 but here is a simple workaround:
Get human reference genome as fasta and suffix all fasta names with _human. Do the same with the mouse genome using _mouse. cat both together and build an index. Then you can later track back whether the alignment was done to human or mouse.
Edit: Be aware though that for those parts of the genomes that have high sequence similarity you will get multimapping reads with MAPQ of 0 (or at least very low).

This is generally the right way to go about it. Different methods will handle problems like sequence similarity or comparing mapping qualities differently, so Nomad420 you should be aware of what different tools do. One place to start looking is the Xenome tool. It is a command line tool that does what ATpoint suggests and classifies reads based on how they map to each reference genome.

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


text:
Heat map of protein expression from normalized abundance
This is a follow up to my question posted here
Processing proteomics data
In relevance to my research, I've been looking for proteomics data (control vs. diabetic) and I found a dataset in the article "Diabetes causes marked inhibition of mitochondrial metabolism in pancreatic β-cells" (figures from here: link).
I would like to know how to generate the following heat map from the normalized abundances reported.
Fold change heatmap:

Abundance:

Figure description:

Heat maps of relative mRNA and protein levels of the indicated genes
in islets isolated from control and 2-week diabetic βV59M mice. Each
box corresponds to a different animal. Colour indicates log2
fold-change. Heat map: TCA cycle. Red, proteins upregulated in
diabetes; blue, proteins downregulated in diabetes. Grey, no change or
not detected. PDK1, pyruvate dehydrogenase kinase 1. Abundance:
Abundance of the indicated proteins, quantified by mass spectrometry,
in islets isolated from control (black, Ctrl, n = 4) and 2-week
diabetic βV59M (white, Diab, n = 4) mice.

I am not able to understand how the log2 fold-change is computed from the abundance.
For instance, if we consider PDK1 (first panel) the abundance measure of control is ~50000 and diabetic is ~200000. Heat map says log2 fold-change (in the figure description added above). What is the fold-change that is referred to? From what is shown, I think this is not diabetic/control fold change; we see two panels (i.e one for control (4 samples) and the other for diabetic (4 samples)) and not a single panel. Since the legend is from (-2 to 2) I think these values are not log2 (abundance) as well.
Could someone please clarify how the heat map has been generated from the abundance measures? Basically, I want to know how to generate the first heatmap figure from abundance plots.

Briefly ... that is absolutely correct - its not clear. The 2 to -2 is log10 (standard log) there is no other explanation. The appendix of the paper used log2 because it stated this (it is weird because natural log is would be the obvious choice in that scenario). This heatmap must be standard log. The authors MUST state the transformation and they do not appear to have done this. I will forward heatmap stuff - as comments.

I do have sympathy with the OPs project - this is an unnecessarily difficult project, which is outside their control.

I got a reply yesterday from the corresponding author on the diabetic vs. control labels of the samples.

Thats good, happy the project is moving forward

I'll delete this comment later - mods complaint - what I'll do put up a python seaborn solution to your data under pandas, but I'm not around until much later in the week. Its  fairly easy to do, so you don't have to hang around.

Hi @M__  Thanks for the comment. I managed to plot the heatmap in python. But I am having some confusion regarding the figure description and the heatmap shown in the article. I'll share a detailed edit today.

The -2 to 2 should be the normalized version of log2 data (maybe mean or median normalization). This is commonly used in heatmaps to show downregulated value vs. upregulated value.


text:
Setting up MaxQuant parameters for processing proteomics data
This is a follow up to my previous question here.
In relevance to my research, I've been looking for proteomics data (control vs. diabetic) and I found a dataset in the article "Diabetes causes marked inhibition of mitochondrial metabolism in pancreatic β-cells" (figures from here: link).
I want to process RAW files from this proteomics experiment. I'm trying to set up the analysis in MaxQuant and I am looking for some help in defining some parameters.
Sample:
PRIDE dataset (https://www.ebi.ac.uk/pride/archive/projects/PXD012979). I want to analyse
8 sample files (please find the file labels and appropriate condition labels below).
Instrument used to generate the data:
Orbitrap Fusion Lumos
Method:
Label free quantification
Modifications:
Groups/condition:
Diabetic vs. Control
Diabetic mice
FL0408_MSQ911_MariaRohm_554_F30
FL0408_MSQ911_MariaRohm_545_F11_544_F33
FL0408_MSQ911_MariaRohm_554_M3
FL0408_MSQ911_MariaRohm_582_M1
Control mice
FL0408_MSQ911_MariaRohm_544_F13
FL0408_MSQ911_MariaRohm_544_F31
FL0408_MSQ911_MariaRohm_582_M13
FL0408_MSQ911_MariaRohm_544_M1
Modifications:
Fixed:
Oxidation(M)
Deamidation NQ
Variable:
Carbamidomethyl(C)
My objective is to get the values of the normalized abundance (kindly find the figure below from the article associated with this dataset, for 4 control vs. 4 diabetic samples) from MaxQuant analysis.
Abundance (shown in the paper):

MaxQuant setup:
I'm not sure what has to be filled in the fields circled in red. Could someone please have a look?
Should the fraction be set to 1? I am not sure what has to be filled for the dataset in the link that I shared above.

Label min. ratio coun -- should this be 1? I am not sure what this means.

Should the following fields be checked?


text:
How can I get the duration of the events in Nanopore?
PacBio has the concept of IPD (InterPulse Duration) which is the time between two detected consecutive sequences in the raw signal. I have been trying to extract this value in the Nanopore, but I cannot find them in the .fast5 files. But this info must be there:
This is from the paper Nanopore sequencing technology, bioinformatics, and applications:
''Each current segment contains multiple measurements, and the corresponding mean, variance, and duration of the current measurements together make up the ‘event’ data.'' Directly quoted from the paper, which suggests we should have the information about the duration in the raw signal.
How can I get the duration, I appreciate it if can anybody help?
Some info I read was confusing to me (e.g. http://simpsonlab.github.io/2015/04/08/eventalign/ which is nevertheless a great source for understanding what is going on), I knew single measurements' speeds are controlled, however, I didn't know the block stride is fixed and it made me so confused at first.
What I found related to my own question:
I guess I was not completely correct in the terminology I was using, which might cause confusion. In nanopore sequencing, translocation duration refers to the time a fragment of DNA/RNA passes through the nanopore.
Also for better visualization, I made this plot after running albacore basecaller:

as @gringer has explained, the block stride value is fixed (=5).
The figure shows a small portion of albacore Event Table:
The x-axis corresponds to the relative time from the start of measurements and the y-axis shows the current level.
Each black dot shows a single measurement from an ONT device; the dots are segmented by a set of block strides (=5) used by albacore to call an event. The red lines show the average of the measurements in one segment. The shaded area in light gray represents the event that resulted in calling a base (one to one). Dark gray areas show events that resulted in
calling two bases (split events). The white area corresponds to events with no called bases (skipped events). The boxes below the x-axis correspond to bases called in shaded areas and the bold base is the base that is called for that event.

I was under the impression that PacBio was one sequence per ZMW. Did you mean two detected consecutive bases?

Unfortunately ONT no longer produce event data as default in their FAST5 output with the most recent basecallers, and they've been quite vague about how to generate it.
The easiest way to fix that is to use an older basecaller. If you tell Guppy v5.1.15 to output FAST5 files, it should include an event table that includes read start time and move information. This might also work with new Guppy, but ONT has told us that the event tables have gone away.
Processing the move tables to get information about signals per base is a bit complicated. Just using the metadata provided in the FAST5 files it's actually impossible, but you can get close enough to match what the ONT basecallers were trained on (which is probably close enough). To get closer than that, you need to process the raw signal yourself (e.g. see here).
Alternatively, you might be able to get that information out by following the tomboy documentation, which produces output that demonstrates signal-level changes:
https://github.com/nanoporetech/tombo
ONT have had various methods for basecalling over the years, and these methods change how the signal is chunked and processed.
When nanopore sequencing was first commercially available, base change points were resolved  down to a single-sample level (as described in Jared Simpson's post which you linked). Or, at least, as accurately as they could at the time. Even with that approach, events were occasionally misclassified such that two consecutive events could represent the same base, or two consecutive events could involve a change of more than one base. To account for this, ONT created the move table, which described how many bases were shifted between one base and the next base. Initially ONT were using pentamers (i.e. 5 base sequences) for the output of their calling models (this had a physical reason, because 5 bases was about as many as would contribute to the instantaneous current), so the move table could have values of 0 [indicating two adjacent events described the same base], 1, 2, 3, 4, or 5. Although 5 was the maximum, because it was a pentamer model, what it actually meant was "there has been a large skip; we have no idea how large."
Unfortunately, sample-level signal parsing is no longer the case, which makes precise event alignment tricky. Most recently, individual signal samples are aggregated into what are called strides, which describe the smallest length of time that the basecaller accepts as representing a single base transition.
For this explanation, I'm going to assume that the directory location in the FAST5 file starts from /baseLoc/, with the FASTQ sequence found in /baseLoc/Fastq. This explanation may become quickly out of date, as ONT change their data and methods frequently.
Raw signal information is found at /baseLoc/../Raw/Signal. The stride length information is found in the block_stride attribute of /baseLoc/Summary/basecall_1d_template. This indicates the number of signal samples allocated for each stride. The move table is found at /baseLoc/Move, and represents the number of base shifts in each stride block (usually either 0 or 1). It sounds like that move table is the statistic you're interested in; multiply the number of 0s before a 1 by the number of samples per stride (plus 1 for the base transition), and that's the base duration in samples, which should be divided by the sampling rate to get the number of seconds per base. Sampling rate can be found in the sampling_rate attribute of /???/channel_id (sorry, don't have a precise location, because I usually just assume 4 kHz... which will probably be changing shortly, from the sound of it).
Some of the signal at the start of the read may have been clipped, and it's necessary to take this into account in order to properly align the move table to the signal. The signal offset is found in the attribute first_sample_template for /baseLoc/../../Segmentation_000/Summary/segmentation, and represents the number of samples (not strides) that should be clipped off the start of the raw signal. There's one more thing that may help for aligning bases to signal, and that's the base offset relative to the move table. The basecaller works in chunks of multiple bases (kmers), so the base corresponding to the move table is actually somewhere in the middle of the called/modeled kmer. ONT doesn't provide any insight into this (it's possible they don't know), so I just try some offsets between 1-5 bases and pick the one that looks the best.
I have a FAST5 annotation R script that may help a little bit more in understanding how the alignment works, see here, specifically the code starting from the comment ## Load data from called fast5 file.
Bear in mind that "number of samples per base" varies in a non-smooth fashion with quite a large dynamic range within an individual read sequence, and could be anything from 1 to over 50 (although I stop believing the basecaller if it reports anything over about 20). It doesn't make much sense calculating it as an aggregate statistic for each sequence, except possibly for detecting really poor basecalling (which is more easily picked up by looking at Q scores in the FASTQ file).

Thank you again! This is by far the best answer and explanation I received. Now the procedure is much clearer to me. After running guppy, now I have the move table, block stride and sampling rate information. Tombo also added some other info, but I think I will move forward with guppy for now.


text:
Convert Mascot mztab file format to XML
Cross-posted here
I would like to convert MASCOT mztab file to XML format. I have been looking for converters and I came across PRIDE Converter but this is no longer available for download from here https://github.com/PRIDE-Toolsuite/pride-converter-2.
Could someone please suggest alternate converters, if available?

There are a number of converters to be found at http://www.pastelbioscience.co.uk/resources/databases.html , although not sure that the one you are looking for is among them. May also be better to ask this question on twitter. There is a large #proteomics community that are very helpful and usually come up with an answer :-)

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.


text:
Tree Building Algorithm that treats gaps as deletions
I'm part of a nanopore sequencing experiment that will sequence several generations of viruses. The intent is to perform directed evolution by putting selective pressure on these viruses and tracking the various mutations that occur. See here for the paper this methodology is based on.
I've been working on a software package that can compile consensus sequences from longread nanopore sequencing (see here). However, there's a wrinkle with generating phylogenetic trees from this type of data. It appears that most tree building software ignores indels. This makes sense, as gaps could easily indicate a lack of coverage from illumina sequencing methods or similar issues rather than betray the existence of an insertion/deletion. However, the nature of our sequencing method is such that a gaps really should be evaluated as mutational changes that we want to keep track of in a phylogenetic tree.
We've explored using NextStrain to generate our trees and visuals, but the tree building algorithms they employ too suffer from this gap assumption (see the issue I raised here). These include IQTree, RAxML and PhyML. I have not had much luck finding an algorithm or software that treats gaps as indels when building a phylogenetic tree. Does anyone know of such an algorithm or package I could employ? I figure I could take a stab at writing my own algorithm, but consider that a Plan B resort. No sense re-inventing the wheel and all that.

Various distance measures can do this given an MSA. See for example https://academic.oup.com/bioinformatics/article/28/4/495/212883. As long as you have a transition matrix that deals with gaps it should work fine. You can build trees with e.g. neighbor joining from distance matrices. References on ML methods for dealing with indels can be found in this paper.

What the OP is saying is correct wft mainstream phylogenetics: the reason is they model point mutations - not alignment reliability. There are numerous methods (> 2 contrasting methods) to do what the OP seeks, its quite easy. There's no need to write their own algorithm and would be discourage (account for lots of parameters). Generic NJ is okish (its pairwise), its not quite what the OP wants. The virus cited in the Cell paper is Sindbis virus (alphavirus) ... very interesting. I think I'm probably going to sit out on this question, but I definitely look forward to seeing the paper. Upvoted.

Speaking personally it will take me a good half day to understand the VEGAS system (Cell paper), which - whilst I get the gist: using viruses as replication vectors to accelerate adaptation, in particularly I think it leverages copy number thresholds of capsid proteins; I don't understand the nuance of the wet-lab system.

I need to do a more thorough testing to determine just how effective this answer may be, but I think I've put together a workflow that is close enough to an appropriate answer that I'd post it here.
In my research I came across gappy, (GitHub here), a piece of software that identifies potentially relevant gaps in a phylogenetic analysis and returns a PHYLIP-formatted file that can be run through most tree-building programs. Gappy accounts for indels, but not substitutions, so the output would need to be combined with another analysis to generate the tree I want.
Fortunately, IQTree allows for a partitioned analysis of mixed data (see here), so by feeding the original alignment AND the gappy output through IQTree you should be able to generate a tree that accounts for both substitutions and indels.
Again, more thorough testing is required, but I tested this possibility using an MSA example taken from this website. Saving the FASTA-formatted example as alignment.fasta, and running that file through gappy2, you get an output file called GAPPY2_alignment.fasta_2-inf.phy. Thus if you create a NEXUS-formatted file like the following:
#nexus
begin sets;
    charset part1 = alignment.fasta: *;
    charset part2 = GAPPY2_alignment.fasta_2-inf.phy: *;
end;

and run it through IQTree using the command iqtree -p example.nex, you get the following tree:

+---ENA|CAA23748|CAA23748.1
|
+-------------ENA|CAA24095|CAA24095.1
|
|              +-------------------------------------------ENA|BAA20512|BAA20512.1
+--------------|
               +-----------------------ENA|CAA28435|CAA28435.1

Now, disclaimer: this is the exact same tree you get when you run alignment.fasta through IQTree, hence why more thorough testing is required. The sample size for this example is quite small, and it's therefore not surprising that accounting for indels does not change the outcome in this scenario. However, I think this serves as a solid baseline workflow, where future updates will be more tweaks than complete overhauls, if that makes sense.
Let me know if you see any glaring flaws in my logic, but I think I'm going to close this question. I may provide edits to this answer in the future if significant tweaks/overhauls wind up being done.


text:
IBM SPSS- 2-years survival of cohort
Over the past few months I have been retrospectively analyzing patient files for my master's thesis in medicine. As the median overall survival of my cohort is low, I have decided that I want to analyze the 2-years survival of my population. However, I have not been able to find any literature on how to perform this analysis on SPSS.
Does anyone here have any experience analyzing the median survival over a specific period of time (i.e., 1-year, 2-years or 5-years) using IBM SPSS?
Thank you in advance!
Leon

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

What you want is Cox regression (proportional hazard) models to perform survival analysis. It is years since I have used SPSS and would encourage consideration for the use of R (Python is a bit more complicated). There is a commercial go-between which is easier to use than R but more computationally efficient than SPSS, but I cannot remember its name at all.
There is a weblink for using Cox regression in SPSS:  https://www.scalestatistics.com/cox-regression.html
I have no idea whether the information is correct, but it looks simple enough to point and click.
Background
A  proportional hazard model is one of numerous survival analysis which model unit risk, e.g. mortality, over time considering the risk e.g. mortality of the baseline population (covariates). Generally a hazard model will comprise a hazard function of the risk of events per unit time against the baseline covariates, such as age, smoking that sort of thing. Thus you wouldn't want to consider the mortality risk of a 25 year old to comparable to that of a 90 year old patient for a given factor, because there's loads of confounders.
At a minimum SPSS could be slow and this is a regression analysis, but you will find out soon enough.

It appears there's a forward prediction being sought here. A simple regression may be able to capture this, but the gist of the question appears to be a precise temporal forward prediction. Precise forward prediction are complicated and whilst there are definitely R packages that do this, at a precision level of forecast e.g. 2 vs 5 year this definitely ain't SPSS, a long way from its capabilities. The OP will need to learn R to do this study, it isn't trivial, or approximate from a hazard style regression. For this style of calculation significant epidemiological and clinical information is required. For example, these are often used in cancer studies, but precisely what the background epi/clinical is will be important as to whether a given R library is relevant.
Summary The answer to the question is SPSS can't do this and learning R is the way forward, because the OP will need to use an established library and this only place I'm aware that such libraries exist. The only alternative is to approximate using classical survival analysis.

Hi there! thank you for your response.
You are absolutely right, it's as easy as that to perform survival analysis. However, this results in a comparison in overall survival, whereas I'm looking specifically for a 2-years survival analysis but I just can not seem to find out how to do this.


text:
Using QCTOOL v2 to process UK Biobank .bgen files - why so slow?
I’m currently using QCTOOL v2 to process imputed .bgen files from UK Biobank, however they seem to be processing very slowly. Is this normal?
My command is pretty basic; I’m filtering out a list of SNPs and samples:
/path_to/qctool \
-g /path_to/ukbXXXX_c22_b0_v3.bgen \
-s /path_to/ukbXXXX_c22_b0_v3.sample \
-og /output_path/ukb_c22_filt.bgen \
-os /output_path/ukb_c22_filt.sample \
-excl-rsids /path_to/snps_rem_c22.txt \
-excl-samples /path_to/samples_to_rem.txt

It is currently processing SNPs at a rate of 1.2/s (e.g. 71169/?,57205.8s,1.2/s).
The computational facility I'm using should not limit the speed of an operation.

Have I made any mistakes in my qctool query?
Or any other ideas on why it's running so slowly?

Alternative suggestions on how to process these files would be appreciated for example, although I would prefer to use QCTOOL perhaps I have to use PLINK.

From the website, https://www.well.ox.ac.uk/~gav/qctool_v2/
"QCTOOL is a command-line utility program for manipulation and quality control of gwas datasets and other genome-wide data"

Hi @Cat could you provide information on computer(s) you are using, e.g. are you using a desktop, laptop or the cloud? What's the RAM of the machine, how full is the hard drive and how big are the files, that sort of thing.

crossposted https://www.biostars.org/p/9519364/

Also whether the files (either input files or output files) are on a local or network drive, that can make a huge difference to the speed.

bgen/qctool2 is just pretty slow in general I've found. For performance, convert to pgen format and use plink2, seems to be orders of magnitude faster as long as you have enough RAM to hold the files in memory,

I contacted the QCTOOL team directly and received this response:

The design of BGEN means that when you subset samples the data has to
be recompressed - this is essentially what makes this slow.  (By
contrast you can subset SNPs very quickly without recompression using
bgenix https://code.enkre.net/bgen.)  It is therefore definitely worth
considering not subsetting samples but using a sample
inclusion/exclusion list instead, if your analysis software supports
that.
If you have to subset and want to use QCTOOL - some things to try are:
I typically use the options -bgen-compression zstd -bgen-bits 8 now,
c.f. https://doi.org/10.1101/308296 this is faster but first check
your downstream software supports zstd compression. Use a map/reduce
type pipeline (i.e. chunk data for re-encoding) - this can be
implemented using bgenix and cat-bgen.

I assume you’re working with the latest version?
Have you tried first stripping the SNPs, then stripping the samples?
Ever since working with the imputed data from the UKB I was never
successful in chopping pieces of bgen files using qctool in a timely
manner. If my memory serves me well, qctool manages to strip SNPs
fairly quickly, but is super slow to remove samples. Try PLINK 1.9 or
2.0.
As for “The computational facility I'm using should not limit the speed of an operation.”,
depending on how disk data moves in/out of slave nodes and how busy
the cluster is, I did saw processes to become I/O starved on large
clusters.
Also, consider that most tools that perform association testing can
utilize SNPs  and samples lists. BOLT, SAIGE, regenie, SNPtest… so
there might not be a need to pre-filter. Put “NA” as phenotypes you’d
like to be “removed” from the testing.

In the end I did not pre-filter SNPs or samples - I set samples to NA within the phenotype file, and used a SNP inclusion list in the second stage of SAIGE with the flag idstoIncludeFile.

When subsetting SNPs (and keeping all samples), I find bgenix to be about 100x faster than qctool.


text:
Subsetting with harmony
Currently merged two Seurat objects together and then ran Harmony for batch correction. Now I want to subset out a cell type of interest do I re run harmony or just do the standard PCA?
scData <- scData %>%
    RunUMAP(reduction = "harmony", dims = 1:14) %>%
    FindNeighbors(reduction = "harmony", dims = 1:14) %>%
    FindClusters(resolution = 0.2) %>%
    identity()

numPC = 14
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.2, dims = 1:numPC)
scData <- RunUMAP(scData, dims = 1:numPC)

Both options give me different results. Not sure what the correct way to it is.

From the Broad Institute website:
https://portals.broadinstitute.org/harmony/
"Harmony is a general-purpose R package with an efficient algorithm for integrating multiple data sets. It is especially useful for large single-cell datasets such as single-cell RNA-seq."


text:
Is there a software that allows to convert a jpeg graph into dataset?
I have a very old graph as a .jpeg picture that I would like to convert into a dataset to be able to operate with it using R/python.
I have no idea if it is even possible.

How would that even be possible? How many datapoints are actually there? In theory, you could extract an infinite number, but the original data only had values for a small subset of all possible points on that line. Given that there are no visible points in your graph and that the quality is extremely bad, any extrapolated dataset would be a wild guess at best.

Looks like Graph Grabber can do the job.
https://www.quintessa.org/software/downloads-and-demos/graph-grabber-2.0.2#download


text:
Error in seaborn plot " Horizontal orientation requires numeric `x` variable"
I am trying to plot a box plot with seaborn with the following code
plot = sns.boxplot(
    y='Min PPI distance',
    x='Synergy_percent',
    color='white',
    orient="h",
    data=synergy_df
)
plot.set(
    xlabel='Synergy (%)',
    ylabel='Min. distance',
);

The data synergy_df looks like
    Combination Min PPI distance    Synergy_percent
0   1003- Camptothecin & 1011- Navitoclax   2.0 22.2
1   1003- Camptothecin & 1017- Olaparib 2.0 6.4
2   1003- Camptothecin & 1022- AZD7762  2.0 40.4
3   1003- Camptothecin & 1025- SB216763 3.0 2.1
4   1003- Camptothecin & 1032- Afatinib 3.0 2.2

After running the above code I am getting the error
TypeError: Horizontal orientation requires numeric `x` variable.
I want to add that Synergy_percent column has many zero values and about 1500 rows; I haven't tried it with a smaller number of rows.
I think the reason for the error is the wide range in the synergy_percent values. The range is between 0 to 86.7. If I use xlim here what range should I use to divide the x-axis in 4 quarters (0-25-50-100)?
Background Seaborn is a Python library for more complex graph plotting than the commonly used matplotlib and is often used for boxplots and heatmaps usually using a pandas dataframe. The official documents are described here.

Very quickly because I'm not around.
I think I know what's going on, but don't have enough time to figure the df. The key thing is seaborn needs to know where columns are x and which are y ... within the data frame . This information appears missing and the OP used y='title', x='title without reference to the dataframe. I'm pretty certain thats the bug here.
import seaborn as sns
import matplotlib.pyplot as plt

plot = sns.boxplot(
    color='white',
    y=synergy_df["Synergy_percent"],
    x=synergy_df["Min PPI distance"]
);
plt.show()

Alternative, which seems more what the OP wants ...
import seaborn as sns
import matplotlib.pyplot as plt

plot = sns.boxplot(
    color='white',
    # use your orientation command here
    x=synergy_df["Synergy_percent"],
    y=synergy_df["Min PPI distance"]
);
plt.show()

Blurb Normally the x-axis is the discrete boundary whilst the y-axis is continuous. In you example the x-axis feels like the continuous variable and you are producing a rotated boxplot - thats perfectly fine and seabourn will do that however its  good to get it working first before rotating the axes.
Final The orientation of the dataframe might possibly be an issue, not sure unless its tried.

If there's a formatting issue the only trick I know is
 plot = sis.boxplot(
 .... # same old
 dodge=True,
 );

I don't know beyond, e.g. imposing percentile ranges onto a boxplot.

Well its good its working. Could you try dodge=True. ?? What I think you want is to tweak the 'box and whiskers' of the boxplot - which I'd have to look up. The other possibility (which I don't think you want) is to make 4 boxplots within the range specified. This is last point is about pandas filtering and would be outside seaborn. I'm pretty certain you just need a simple tweak within boxplot method. However, beyond the dodge trick I dunno, but dodge does clean stuff up.

Good to know; I guess that means data=X doesn't have the same effect in seaborn as it does in R.

I'm not familiar with seaborn, but I'd guess from the error that there's at least one text value inside the Synergy_percent field that is encouraging the rest of the values to text. Either that, or spaces in the input file are messing up the parsing of the columns.


text:
How to do enrichment analysis with scipy.stats.hypergeom.sf function
I want to do pathway enrichment analysis in python using scipy.stats.hypergeom.sf function to select pathways that are significantly enriched in at least one tissue.
my data:
Enrichment  CANCER_TYPE PATHWAY
11.6    Breast Carcinoma    Cell cycle/DNA replication
2   Breast Carcinoma    Protein stability and degradation/DNA replication
2.8 Breast Carcinoma    EGFR signaling/DNA replication
4   Breast Carcinoma    p53 pathway/DNA replication
0   Breast Carcinoma    PI3K/MTOR signaling/DNA replication
1.2 Pancreatic Carcinoma    Mitosis/RTK signaling
3.3 Pancreatic Carcinoma    EGFR signaling/RTK signaling

what is "Enrichment"? usually this kind of data comes in the form of a list of genes with annotations and some sort of background set of annotations for genes in a genome. At a minimum, we would expect some count of genes/gene annotations, which is necessary to get the numbers to plug into a hypergeometric test. Different tools will count annotations differently, perhaps tell us more about where your data table comes from?

Could you please add some more detail to your question? Context helps a lot. It's great that you've included input data, but we don't know how that data is generated, or what the 'Enrichment' field means.


text:
How can I call one file from a channel that's composed of several files in Nextflow?
In Nextflow I have a process that creates a channel that holds several files in a tuple. I'm trying to call only one of those files as input for a subsequent process.
So basically if I run channel_name.view() the channel shows [/home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bestTree.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bipartitions.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bipartitionsBranchLabels.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_bootstrap.tree, /home/nextflow/work/e4/bc47bceb4569cffc01bd76e49245e9/RAxML_info.tree] . I need just one file so I'm using the full name in the input section of the process.
input:
file 'RAxML_bipartitionsBranchLabels.tree' from channel_name

That does not find the file. On the other hand, if I use .flatten(), like this:
input:
file 'RAxML_bipartitionsBranchLabels.tree' from raxml_out_ch.flatten()

It works but reads some files before reaching to the one that matches and so the process results in error even after completing the task for this file properly and gives several outputs (for each file it looks into), some of them empty and one with the result.
[7c/2971f1] Submitted process > tree_visualisation (4)
[8d/dd129b] Submitted process > tree_visualisation (1)
[cb/0f91f5] Submitted process > tree_visualisation (5)
[6c/fb23c8] Submitted process > tree_visualisation (3)
ERROR ~ Error executing process > 'tree_visualisation (5)'

Caused by:
  Process `tree_visualisation (5)` terminated with an error exit status (1)

Command executed:

  Rscript tree.R RAxML_bipartitionsBranchLabels.tree phylogram.png

I'm looking for a way to filter this input so that the process only runs for this file and doesn't try with the others.
How can I do that?

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.

Be more specific, is there a certain text pattern that this file matches? Is it always the first or the last file? Is it a file with an extension? Does the channel contain 1 tuple with multiple files or multiple file entries as separate emissions?

That is an unusual setup, and it will likely end up filtering the channel for that specific file. Please provide a minimal reproducible example.

If, in your input definition, you specify a string to represent your filename and provide multiple file inputs as a collection, a numerical suffix is appended to each filename:

When a target file name is defined in the input parameter and a
collection of files is received by the process, the file name will be
appended by a numerical suffix representing its ordinal position in
the list.

When you use the flatten operator, the collection is flattened such that each entry becomes a separate entry in the output channel. So for each of the five 'tree' files, five tasks are spawned: one for each 'tree' file. These files are each localized using the filename specified: 'RAxML_bipartitionsBranchLabels.tree`. This usually isn't what you want to do.
What you're looking for, I think, is some way to select the file called 'RAxML_bipartitionsBranchLabels.tree' from the initial collection of files. There are many ways to do this. For example if you wanted to flatten the output channel, it's possible to select the file you want using the first operator and a closure:
bipartitionsBranchLabels = raxml_out_ch.flatten().first {
    it.name == 'RAxML_bipartitionsBranchLabels.tree'
}

However you should only seldom need to do this. If your pipeline requires just one of these files in a downstream process, my suggestion would be to alter the upstream output declaration to include a separate output channel for this purpose:
process RAxML {

    output:
    path 'RAxML_*.tree' into raxml_out_ch
    path 'RAxML_bipartitionsBranchLabels.tree' into raxml_bipartitionsBranchLabels

    """
    touch RAxML_bestTree.tree
    touch RAxML_bipartitions.tree
    touch RAxML_bipartitionsBranchLabels.tree
    touch RAxML_bootstrap.tree
    touch RAxML_info.tree
    """
}

Thank you Steve. The last option you propose wouldn't really work with that program and I tried the others without success. What I did was to use .collect() in the input of the following process, this way it turns the outputs (of one channel) into an array (or tuple) and then I used def = ... within the script seccion to assign each file in the array to a letter.

No worries @joselopez. The last option should be tool agnostic unless I'm mistaken. It's difficult to say what the right solution is here without the details. Calling collect() seems redundant, since 'channel_name' (i.e. 'raxml_out_ch') already provides an array (java.util.ArrayList). I think you could also do a def branch_labels = yourfilelist.find { it.baseName == 'RAxML_bipartitionsBranchLabels' } in your script block to get the file you want. Note that this solution will localize all the files in the array, which may not be what you want. Your question mentioned you only needed the one.


text:
Parsing pre-2007 SMILES string
How would one parse the SMILES string
 BrC[2]:C[3]:C(:CH:CH:CH:@2):CH:CH:CH:CH:@3

I rely on tools like rdkit and OpenBabel to parse SMILES, but both tools aren't able to parse this string.
More specifically, this SMILES string comes from the supplementary information of a 2007 paper on solubility prediction. After doing some searching, it seems that this paper was received prior to the 2007 standardization of OpenSMILES (this might be why canonical parsers like rdkit have trouble with it).
After searching a bit more, I found an older PubChem sketcher that was able to convert the SMILES in question into a more recognizable format
 (BrC1=C(C(=C(C2=C1C(=C(C(=C2[H])[H])[H])[H])[H])[H])[H])

(parsing this was successful via rdkit, see molecule image below) but still need a way to convert similarly-formatted SMILES strings into the OpenSMILES standard/canonical format en mass.

Are the brackets dynamic? Thus one bracket pattern for one structure a different one for  a different structure or are they always as stated in the second expression? Hope that makes sense - it does make a difference to the answer.

I don't think they're dynamic but am not very familiar with the SMILES format specification.

"Welcome to Jurassic park" (theme music starts playing). Wow, that is proper ancient. And more similar to a SMARTS than a SMILES!
Your outputted version is kekulised and a machine-made mess —I don't think that's a good option. It ought to be something like [Br]c1c2c(ccc1)cccc2 if keeping your input order, but I think the canonical order is : [Br]c1ccc2ccccc2c1 (unchecked, simply typed out while counting carbons aloud). As you can see the standard isn't too dissimilar, so one could write a conversion schema on a push —although I am not sure if it's worth it.

@MatteoFerla Do you know of any software tools that can execute this conversion? If not, how would you go about converting the older format to a more standard one?

I have not encountered pre-2007 SMILES nor do I know of any tool to do so. Few sites would still exist 15 years later. C++ pre-11 generally compiles —pre-98 is the nightmare. Scripting modules may exists, but in Perl, not Python... So CPAN would be the place to search. But the question is whether this is worth it... after all your example molecule is bromonaphthalene.

I see, thanks. My problem is not parsing this particular molecule but rather a whole host of similarly-formatted ones (around 2k), so an automated parser would definitely help!

I assumed as much, I meant it is a rather small compound and there are hordes of catalogues. However, I had a look at your source and it has a name field. Why not just search the name with PubChem PUG API to get the SMILES?
PS. If you do find a Perl module in CPAN that does do the conversion or a similar converter do post the answer!

I couldn't find a conversion module, but the name field trick worked - thanks! Am posting more details as an answer.

I couldn't find a way to directly convert these SMILES strings to a more canonical format. However, thanks to @Matteo Ferla, I was able to get around the problem entirely by looking up the chemical names  and resolving the canonical SMILES programmatically.
More specifically, I used the CAS no column given in the supplementary information of the referenced paper to pull the SMILES via Chemical Identifier Resolver. See this answer for an explanation on how to do so.
EDIT: Some of the CAS nos couldn't be resolved, so I manually converted them using the PubChem sketcher tool mentioned in the question. There are quite a few (~160) of such molecules, so if anyone needs the full parsed dataset but doesn't want to spend the time manually parsing those molecules, let me know.

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.


text:
different results of RMSD for the same trajectory file
I want to analyze the RMSD of my trajectory files using MdTraj.
I used two different codes for this purpose:
the first code:
rep = 1
traj_name = "step7_%d.dcd" %rep
rmsd_total = []

while (os.path.isfile(traj_name)):
    t = md.load(traj_name, top = 'step5_input.pdb')
    top = t.topology
    protein_sel = top.select('backbone')
    protein_traj = t.atom_slice(protein_sel)
    rmsd = md.rmsd(protein_traj, protein_traj, 0)
    rmsd_without_zero = [i for i in rmsd if i != 0.0]
    rmsd_total.extend(rmsd_without_zero)

    rep+=1
    traj_name = "step7_%d.dcd" %rep

and the second code:
traj_1 = md.load('step7_1.dcd', top='step5_input.pdb')
traj_2 = md.load('step7_2.dcd', top='step5_input.pdb')
traj_3 = md.load('step7_3.dcd', top='step5_input.pdb')
traj_4 = md.load('step7_4.dcd', top='step5_input.pdb')
traj_5 = md.load('step7_5.dcd', top='step5_input.pdb')
traj_6 = md.load('step7_6.dcd', top='step5_input.pdb')
traj_7 = md.load('step7_7.dcd', top='step5_input.pdb')
traj_8 = md.load('step7_8.dcd', top='step5_input.pdb')
traj_9 = md.load('step7_9.dcd', top='step5_input.pdb')
traj_10 = md.load('step7_10.dcd', top='step5_input.pdb')
traj_11 = md.load('step7_11.dcd', top='step5_input.pdb')
traj_12 = md.load('step7_12.dcd', top='step5_input.pdb')
traj_13 = md.load('step7_13.dcd', top='step5_input.pdb')
traj_14 = md.load('step7_14.dcd', top='step5_input.pdb')
traj_15 = md.load('step7_15.dcd', top='step5_input.pdb')
traj_16 = md.load('step7_16.dcd', top='step5_input.pdb')
traj_17 = md.load('step7_17.dcd', top='step5_input.pdb')
traj_18 = md.load('step7_18.dcd', top = 'step5_input.pdb')
traj_19 = md.load('step7_19.dcd', top = 'step5_input.pdb')
traj_20 = md.load('step7_20.dcd', top = 'step5_input.pdb')

combine = traj_1 + traj_2 + traj_3 + traj_4 + traj_5 + traj_6 + traj_7 + traj_8 + traj_9 + traj_10 + traj_11 + traj_12 + traj_13 + traj_14 + traj_15 + traj_16 + traj_17 + traj_18 + traj_19 + traj_20

protein_sel = combine.topology.select('backbone')
protein_traj = combine.atom_slice(protein_sel)
rmsd = md.rmsd(protein_traj, protein_traj, 0)

I believe these two codes should result in the same RMSD. However, I get two different RMSD diagrams out of these two codes?
Can anyone help me what is the problems with my first code? (since I think the second code is more accurate while is more naive)

Hi @amirre this has been authoritatively answered please kindly consider accepting the post and helpful for the other site members.

It's a reference issue. The second term.
Your first block makes a list of numpy arrays which are the RMSD against frame  0 of each of the stacks, while in your second you make a concatenation of frames, so your reference frame 0 is common.
The first code block could be smartened by splitting (abstracting) the file valid name enumeration and the reading/analysis. Namely, generate a list of valid filenames and then analyse. This separation would allow you recycle more easily the code: in the seocnd code block you have the file list problem again, which is instead resolved by a copypaste job —very inflexible. A one liner list comprehension (or, less so, a functools.reduce) is a nice option.


text:
Flux balance Analysis, Result interpretation, solver Solver Heuristics
i am new to modelling of cellular networks, and I have two basic questions, regarding the result interpration.
I am using the COBRAToolbox to do a flux balance analysis of a whole genome scale model of a bacteria. To get things started i replicated the results published in a paper, I have even been in correspondence with the author of the paper which happens to be the person who created the model i am using. So I am 100% sure that all the boundary conditions exchange fluxes and so on are as in the publication. I was able to reproduce this results and they match the experimental data.
The published results specifically mention the flux of a metabolite X of interest and compare the reaction rate of X to literature values and to laboratory experiments.
I did a FVA of X given all the boundaries, the soultion space is between 0 - 15 [mmol/g_dw/h], and the result when doing the FBA is aorund 7 [mmol/g_dw/h]?
(I tried out multiple solvers, showing the same variation)
My questions:

I understood the FVA as a analysis tool which shows me the possible solution space for a given flux , given a system with all limitations and boundaries.

So, how is the solver deciding which of the possible values are displayd, which heuristics are used?

If the scattering for the possible solution space is so huge, what is the quality measure for this results?

Thanks in advance.
Best Regards

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.


text:
Sliding window optimisation
I would like to check transposable elements through a genome with a sliding window approach but I am not sure about which window and step size to use.
Is there a way to optimise the window size and step size when doing sliding window based operations? I mean how to choose which sizes to pick?
Thanks!

This question appears sparse on detail. What genome is this? What database of genomes? Traditionally a lot of work on transposable elements was performed on Drosophila. Also what is the algorithm for identifying transposons, e.g. via LTRs, blasting against a transposon database, or looking for the transposase protein or something else?

This is rabbit genome (Orycun2.0). And basically I would like to check the correlation between transposable elements and nuclear insertions with mitochondrial origin (numts). In the corresponding literature there are informations about the window sizes used but not about the step sizes.

I will reiterate @M__'s point: the sliding window technique refers not to a single algorithm but to a family of algorithms and procedures. It is very difficult for us to give you useful feedback if you don't provide more details on the tool(s) or approach(es) you are using.

I just would like to repeat the same experiment as they did in this great article:
https://www.mdpi.com/2075-4450/11/10/680/htm# in the 4th figure. Where the authors correlated the density distribution of numts and TEs. They provided the window size but not the step size.


text:
Minimap2 detects too many 16S sequences in metagenome
I'd like to extract bacterial sequences from a metagenome where ~99% of contigs are from a host insect. My current protocol minimaps every 16S sequences present in the sample, then blasts all those sequences to create a library of full genomes that I minimap against. I'm using the full metagenome as my reference file. I wasn't sure if that makes a difference, but figured I'd mention it. I ran the following command to retrieve 16S sequences that are present in the sample:
minimap2 All_reads.fastq 16S_reference_library.fasta > 16S_map.paf

Minimap2 retrieves every entry in 16S_reference_library.fasta. It's highly improbable some of those sequences are present, or that 25,000 species are present at all. Are there any parameters I could adjust in Minimap2 that might help with this?
Alternatively, I'm trying to implement Krakenuniq and Centrifuge, but I'm concerned that, because I don't have a reference genome for the host insect, they'll have a very high false positive rate.

Do you mean to be using your reads (fastq file) as the reference rather than the query? Normal syntax is minimap2 reference.fasta query.fastq. What do you mean by "retrieves"? How good are the alignments, are they *? Can you show us a snippet of the PAF?

16S is a very conserved sequence, which is why it's used for targeted phylogenetic analysis; it makes it easy to amplify and analyse. Unfortunately that conservation is an issue with minimap2, which is built around the idea that matching scattered subsequences within a sequence is good enough for identifying matches to the entire sequence. With a highly-conserved sequence, there's a higher chance that subsequences will be identical, meaning that diverged sequences will appear to be the identical from the perspective of the mapper.
Minimap2 can be tweaked to make it less sensitive and more specific, but you'd be better off using a different mapping tool, especially one designed for phylogenetic analysis and conserved sequences. I recommend kraken2.
pre-built 16S indexes for kraken2 can be found here:
https://benlangmead.github.io/aws-indexes/k2
If you really want to map specifically to your own generated metagenome, then you could try LAST (with results filtered through last-map-probs), which I have found to be a little bit better in dealing with specificity issues.
Aside: I notice that you've stated "bacterial sequences from a metagenome". If these are, in fact, metagenomic reads (and not targeted 16S reads), then you would be better off using one of the larger metagenome bacterial databases  (e.g. plusPF or plusPF-8, depending on computer memory) instead of concentrating only on 16S-matching reads. This will provide much better resolution; 16S only really works down to the family level.

You are taking reads as the reference. Given the similarity of 16S, you will find hits for most 16S in other species. You have to filter out poor hits. This is the faster way if you don't have many contigs.
Alternatively, switch the order of reference and reads. But in this case, you need to use a high k-mer occurrence cut off as most of your k-mers will occur >10,000 times. Something like:
minimap2 -c -f30000 -N100 ref-16S.fa reads.fq

This will be fairly slow as there are too many hits in the database but it will be easier to process the output.
You can try both and see what works for you.


text:
What is & how to solve File error: my.xml.state (Remote I/O error)?
I caught the next exception during my phylogeographical analysis in BEAST 2 with GEO_SPHERE. What could be the reason? & how to evade this in the future?
      ...
      856000000     -3662.2647      5969.5577     -9631.8225 42m16s/Msamples
      857000000     -3522.3540      6105.9849     -9628.3389 42m5s/Msamples
      858000000     -3621.7922      6022.9359     -9644.7282 41m54s/Msamples
      859000000     -3463.1050      6107.9766     -9571.0817 41m43s/Msamples
      860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples
File error: VP1_test.xml.state (Remote I/O error)
java.lang.RuntimeException: An error was encounted. Terminating BEAST
        at beast.app.util.ErrorLogHandler.publish(Unknown Source)
        at java.logging/java.util.logging.Logger.log(Logger.java:979)
        at java.logging/java.util.logging.Logger.doLog(Logger.java:1006)
        at java.logging/java.util.logging.Logger.log(Logger.java:1029)
        at java.logging/java.util.logging.Logger.severe(Logger.java:1776)
        at beast.app.beastapp.BeastMain.<init>(Unknown Source)
        at beast.app.beastapp.BeastMain.main(Unknown Source)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.base/java.lang.reflect.Method.invoke(Method.java:566)
        at beast.app.beastapp.BeastLauncher.run(Unknown Source)
        at beast.app.beastapp.BeastLauncher.main(Unknown Source)

Unfortunately, googling this error with "BEAST" gave me nothing. If without it, there are a lot of irrelevant problems.

Update
I have just found out that this error happened simultaneously onto all my nodes (different processors, RAM, & common hards & SSDs). All logs of BEAST 2 analyses have it. What is stranger, some of the analyses pressed on after this & someones were broken. To the present I got more interested not in the broken ones but could this exception spoil those which continued their work?

You can restart this at the place you left 'cause is Beast2.
  860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples

However, this is almost 1 billion replications and should be sufficient for convergence. The fact it ends on a round number might infer a algorithm limit is in place. If the Beast2 was restarted from this point and it abruptly ended this would be a very likely reason.
The alternative is an error external to the algorithm such as a RAM error - it has been running a long time.
Whatever the reason I would simply perform ESS statistics and convergence diagnostics on what you have now. There's no real reason to continue if it has converged for 10e6 replications maximum and this is much, much bigger than that.

The terminated processes have not spoiled the work again you can simply input the terminated chains into the restart for Beast2. This is because Beast2 uses checkpoints. Brief information on how to resume from the last checkpoint is highlighted here https://www.beast2.org/2014/04/07/checkpointing-tricks.html. Convergence checks are important.
The situation is only is only catastrophic for Beast1 where you would have needed to start from the beginning.

Unfortunately, big number of parameters demands big number of iterations. I have several thousands of leaves, some sequences are the same but with different dates & locations. Moreover, continuous phylogeography with complex region shapes is really a hard task.


text:
How to convert Foldseek .m8 alignment files to FASTA format?
I want to use sequence alignments found by Foldseek in PROSS. I know there are other ways to get alignments (also tried with results of HHPred), but I had in mind that Foldseek may find remote homologs that would still be structurally relevant.
Currently Foldseek only gives alignment results in m8 format, however PROSS only accepts FASTA alignment files.
I tried looking into BioPython but it seems it doesn't open m8 files, only m9. I also tried importing the alignment in Clustal X but it can't read it.
I could probably use Foldseek CURL API to get the alignment as json and make a script to write the results in FASTA format, but this seems like a heavy handed solution.
What would be other tools I could try to convert m8 to FASTA?

tl;dr
Use an awk script (below) or "one-liner"
awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta

Caveat emptor: I've tested this on my own Foldseek outputs and it works for my pedestrian purposes of multisequence alignment with COBALT. I haven't tested it with PROSS but expect it should work for you too.
Long answer:
Using this definition of m8 (which is a modification of the BLAST m8 format):
https://github.com/steineggerlab/foldseek#output-search
And this definition of FASTA:
https://en.wikipedia.org/wiki/FASTA_format
We see we want to take lines like this (I've added [...] here to show what I've cut off for brevity):
job.pdb_A   Some matched structure name 51.100  311 138 8   2   305 25  328 1.000   2.790E-35   1272    306 334 VDVSVSIFINKIYGVNTLEQTYKVDGYIVAQWTGKPR--KTPGDKPLIVENTQIERWINNGLWVPALEFINVVGSPDTGNKRLMLFPDGRVIYNARFLGSFSNDMDFRLFPFDRQQFVLELEPFSYNNQQLRF--SDIQVYTENIDN[...] -27.488,-8.287,-35.993,-26.838,-4.500,-36.421,-25.278,-3.958,-32.938,-22.398,-6.436,-33.555,-21.375,-4.861,-36.892,-21.373,-1.326,-35.385,-19.237,-2.585,--18.044,-1.289,8.073,-14.345,-2.072,7.367,-15.099,-5.282,5.328,-16.872,-3.358,2.522,-16.876,-5.018,-0.77 [...] MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]   1918945 Vibrio spartinae

and map them to entries like:
>some matched structure name
MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]

We can do this with an awk script (The awk manual is your friend here https://www.gnu.org/software/gawk/manual/html_node/Long.html). Put the following in a file somewhere handy and call it using awk -f m8_to_fasta.awk my_m8_file.m8 > my_fasta_file.fasta
you may want to combine this with head, tail, and/or grep to reduce the number of sequences that end up in your output.
#! /bin/awk
{
        name = "> "
        # start on the second field, foldseek will put the job name in $1
        field = 2
        # Check each field until we see numbers
        while ( $field + 0 != $field ) {
           name = name " " $field
           field = field+1
        }
        # name now contains the name string which is the first line of a fasta entry

        field = NF
        # Now find the matched sequence
        # Sometimes an entry ENDS in a number, so we strip that
        if ( $field + 0 == $field ) { field = field -1 }
        while ( $field + 0 != $field ) {
            field = field - 1
        }
        # field is now pointing at some undocumented number. Decrement once more to get the sequence
        field = field - 1
        # field is now pointing at the sequence
        seq = $field

        #Now I/we want to sanitize the sequence to remove any gaps
        gsub(/-/, "", seq)

        print name
        print seq
    }

Or we can coerce this into a monstrous one-liner:
 awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta


text:
Missing value imputation method for gene expression data
I am new to working with gene expression data sets and am wondering what is the most standard or best way to impute missing values in a gene expression data? I got mine from the GEO database and the accession number is GSE70947. This is a breast cancer data set with a total of 296 samples (148 normal/148 cancer) and there are 62976 features (genes).
The original data is only log2 normalized by the experimenter.
When I first queried the data from GEO database using GEOquery from Bioconductor, the data set had NAs in it and a lot of them -- I found that there are 729494 of them in the data set.

sum(is.na(gene.expr))
[1] 729494

Now, I am looking online to see what are the most common imputation methods for gene expression data, but there are so many (local, global, hybrid, knowledge assisted methods) and I don't really know which is more suitable for my data set. I am just wondering what is the best way to impute the NAs. Should I just impute the NAs to 0s?

Once I imputed NAs to 0s and that was a mistake. NA is 'no available' record and it can be anything. So you can work with logs, do na.omit of your z-scores or download raw data and process it by yourself with bioconductor. It's Agilent-028004 SurePrint G3 Human GE 8x60K Microarray, so it should be a package to deal with it. Then you can do anything with normalized signals.

What about Bayesian approaches like mc-stan.org? Those generally provide good approaches for working with missing data.

@Adamm  You said, "It's Agilent-028004 SurePrint G3 Human GE 8x60K Microarray, so it should be a package to deal with it". How so? Where can I find it?


text:
What is the different between the I-TASSER, phyre2, SWISS-model in the 3D tertiary structure?
What is the difference between the I-TASSER, phyre2, and SWISS-model in the 3D tertiary structure?
How do they get the results? When I did a prediction, I got a similar result for the highest template in I-TASSER and Swiss, while it was different with Phyre2. What does that mean?

This is a common Q.
Briefly, you have some techniques here that get used separately or combined:

Ab initio: assemble the protein from scratch, often via a library fragmented peptides (3-mer, 6-mers, 9-mers) based on real structures.
Template threading: given a protein structure alter it to make it like a give sequence. SWISS-MODEL, Rosetta's SimpleThreadingMover, ModBase, Phyre template threading etc. This is why some of these keep ligands and deal well with multimers. For missing loops in the template these are filled by ab initio methods.
Hybridisation: give a bunch of models mix and match them for the best energy. I-Tasser, RosettaCM (which uses SimpleThreadingMover), Phyre (which uses Phyre template threading)
Ab initio with contraints, such as evolutionary covariance. EVFold, AlphaFold2, trRosetta. Often mixed with template based approaches.

Wikipedia has more, but papers are very helpful despite the car-salesmanship.
One thing to note is that protein can adopt different conformations. The AlphaFold2/ColabFold notebooks give multiple models: these may be legit alt conformations, e.g. inactive vs active etc etc
Here is a SwissModel vs AF2:
Different structures of same protein from two softwares
Other helpful As:

How to make 3d model of a protein that not exist in PDB?

Phyre2 vs ITasser, completely different models generated

How reliable is Alphafold2 for folding random-coil peptide —discussion of technical details going wrong... Personally, I had Phyre2 made an alignment mistake which altered the protein


text:
Reference variant detected as altered one in bam file
I received (from manufacturer) several .bam files and I used four callers (samtools, freebayes, haplotypecaller, deepvariant) to find some sequence variants. In obtained .vcf files, I took a closer look to some calls. I found interesting, homozygous one rs477033 (C/G Ref/Alt) with flag 'COMMON=0' and very low MAF. I also checked the .bam file (in IGV) and this rs itself in ensembl database. It looks like this:

The problem is: in this case ref and alt are switched. The G is reference (ancestral) allele and C is altered. Nevertheless in the .bam file it's interpreted in opposite way. Of course it's false positive variant - there's no mutation there at all. All variant callers returned it as high quality variant. My. bam files are full of them.
What went wrong in this case? Problem with mapping?
How can I fix this thing, to get only novel and rare true positive variants?

[link] https://drive.google.com/file/d/1eheG5MqS_a1xK5HSqsNJ0Tito_WdL_YV/view?usp=sharing

A reference sequence can obviously not represent all populations,  as shown above in the linked chapter (• page 30 •),  a reference allele may be the minor allele in the population you are studying ( in some cases )

The ancestral allele is not always the reference allele. Ref and alt aren't switched, ref is just rare in at least the population you're looking at (it appears to be rare in humans as a whole). So, the C is indeed reference and G is indeed altered, it's just that people homozygous for the G are the overwhelming majority. This is evidenced by the rs477033 variant being named NC_000001.11:g.165890393C>G. So, in short, nothing went wrong, all tools functioned correctly.
You'll notice a bunch of these over time, either make a database of them (or find one, someone may have made one) and filter them out or see if maybe the 1000 genomes reference genome has some of these modified (or maybe the T2T genome that was recently published).

I think that I found the solution and I'd like to add it here. It might be helpful for others. I just simply performed additional annotation of my .vcf files using 1KGenomes data. I've found it here: here Then each variant is described in INFO by AA which means ancestral allele. If this allele is different comparing to my REF I filter it out.


text:
Biopython reads my tree eternally long
I have a nexus tree (1332 taxa) with a lot of additional data. When I tried to read it through tree = Phylo.read(treepath, "nexus"), my kernel got eternally loaded.
If I abort the process, I get the following message:
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Input In [95], in <cell line: 1>()
----> 1 tree = Phylo.read(r"D:\Downloads_hard_disk\Job\NAs_Hepatitis_A\Results\experiments\test_right\test_right_ExpPop_3\test_right_ExpPop_3.biopy.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1181, in Nexus._tree(self, options)
   1179     elif special == "W":
   1180         weight = float(value)
-> 1181 tree = Tree(name=name, weight=weight, rooted=rooted, tree=opts.rest().strip())
   1182 # if there's an active translation table, translate
   1183 if self.translate:

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:82, in Tree.__init__(self, tree, weight, rooted, name, data, values_are_support, max_support)
     80 # there's discrepancy whether newick allows semicolons et the end
     81 tree = tree.rstrip(";")
---> 82 subtree_info, base_info = self._parse(tree)
     83 root.data = self._add_nodedata(root.data, [[], base_info])
     84 self._add_subtree(parent_id=root.id, tree=subtree_info)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Trees.py:123, in Tree._parse(self, tree)
    121 elif not incomment and tree[p] == ")":
    122     plevel -= 1
--> 123 elif tree[p:].startswith(NODECOMMENT_START):
    124     incomment = True
    125 elif incomment and tree[p] == NODECOMMENT_END:

KeyboardInterrupt:

I am interested what is the specific reason of this: either the ineffective algorithm or a bug in Biopython? Is there a way to read the tree as Bio.Phylo.Newick.Tree object? Or should I completely use another package? In the last case, could you recommend effective one?

The reason is the RAM has been exceeded for Biopython's capabilities and you are encountering a 'RAM-bottleneck'. I've listed the options from 1 or 3 and recommend option 3c.
The options are use:

machine with much greater RAM
acceleration coding techniques in Python (such parallelisation)
Beast tools written in Java

I assume this calculation is being performed on your laptop because in theory <1500 isn't much, but agreed there are ~1500 taxa.
Option 1 I personally don't think you are very far into the bottle-neck. Certainly a Cloud computer (e.g. GCP) could easily perform this calculation and maybe your the headnode of your local cluster. Usually the headnode has the most RAM.
Option 2 Code acceleration is worth keeping in mind because this achieves >> 1 log performance. Parallisation will solve this using an optimal chunk-size. However, large amounts of code can be needed. Code optimisation - Python is very sensitive to this because it uses C-libraries - for example, use map and reduce and avoid ANY regex in your Python code because IMO its a key weakness of Python. Pandas iterrows is famously slow. If the Biopython library is the problem code efficiency isn't an option (unless you want to rewrite Biopython). Label compression or header recoding techniques work, I'm not sure whether this has been done in Beast output already(?).
Option 3 The tools for dealing with this which are written and optimised to achieve rapid output are:

Tracer - does lots of analysis on MCMC output. Its really geared towards assessing convergence, but it can do a lot of stuff.
TreeAnnotator - makes a summary tree from the MCMC output you have

TreeStat2 <- I think this is what you are looking for.


text:
Protein-Ligand in PDB
I am a newbie in bioinformatics. In PDB, while searching for proteins, I sometimes get protein with ligand, such as PDB ID 1J1Z. In cases I get pure protein structure, in some cases I get protein-with ligand.
I am having hard time understanding why I can have ligand in PDB structure. I am PDB is supposed to give me pure protein 3D structure. Then how can we have ligand in protein structure??
And sorry if my question is stupid.

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

In some protein crystallography studies, the protein crystals are grown with the substrate/ligand/inhibitor added, thus the resulting structure contains both the ligand and the protein atom coordinates. This enables studying the protein ligand complex to understand binding, specificity, mechanism of catalysis...
On the other hand in some experiments protein crystals are grown without the ligand, thus the resulting structure contains coordinates just for the protein atoms.
On the third hand, the lucky one, your protein, overexpressed in its natural host [sometimes even in different one], co-purify with its own natural ligand, so the resulting structure contains coordinates for both the protein atoms and the sought for ligand.
To understand ligand binding the protein structure without the ligand should be compared to the structure of the protein-ligand complex.
Many published papers report both types of crystals.
For example the crystal structure of pectin methylesterase in complex with hexasaccharide VI: 2NTP is from the publication https://doi.org/10.1038/sj.emboj.7601816 which is the source of six other structures. These structures include the enzyme with different substrates, as well as the active site mutants (D178A) all of which were prepared to understand the molecular basis of substrate specificity. The conclusion after studying all of the mentioned structures is nicely summed in the sentence from the abstract:

The preferential binding of methylated sugar residues upstream of the catalytic site, and demethylated residues downstream, drives the enzyme along the pectin molecule and accounts for the sequential pattern of demethylation produced by both bacterial and plant PMEs.

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:
Why must a maximal non-branching path be a contig?
The following is from Bioinformatics Algorithms:

Fortunately, we can derive contigs from the de Bruijn graph. A path in a graph is called non-branching if in(v) = out(v) = 1 for each intermediate node v of this path, i.e., for each node except possibly the starting and ending node of a path. A maximal non-branching path is a non-branching path that cannot be extended into a longer non-branching path. We are interested in these paths because the strings of nucleotides that they spell out must be present in any assembly with a given k-mer composition. For this reason, contigs correspond to strings spelled by maximal non-branching paths in the de Bruijn graph. For example, the de Bruijn graph below, constructed for the 3-mer composition of TAATGCCATGGGATGTT, has nine maximal non-branching paths that spell out the contigs TAAT, TGTT, TGCCAT, ATG, ATG, ATG, TGG, GGG, and GGAT. In practice, biologists have no choice but to break genomes into contigs, even in the case of perfect coverage (like in the figure below), since repeats prevent them from being able to infer a unique Eulerian path.

But this quote doesn't make sense to me:

We are interested in these paths because the strings of nucleotides that they spell out must be present in any assembly with a given k-mer composition.

Is it not possible that we simply have a missing read that would have otherwise branched this path? How do we know that a maximal non-branching path must be a contig?

The definition of contig is derived from assembly algorithms. A contig is not necessarily present in the genome. PS: i.e. "contig" is algorithmically defined, not biologically defined.

Is it not possible that we simply have a missing read that would have otherwise branched this path?

Yes, it is possible. Then you are likely to get a misassembled contig.

By the way, a maximal non-branching path in the assembly graph is more often called a unitig.


text:
Can I use a benchmark to choose a clustering or module detection method for a gene dataset?
The paper A comprehensive evaluation of module detection methods for gene expression data says that: "We first want to provide an overview of the characteristics and performance of current module detection methods to guide the biologist in their choice".
My problem is that I have no gold standards for my gene dataset and I would like to know if I can use the paper results to choose a module detection method without the need of external validation measures and how can I justificate this generalization? I am not finding the justificative about how can I use this guiding.

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.

How does your dataset compare to the dataset in the paper? Without that information we can't answer the question (as noted also in your previously posted SE bio question).