# Question
Title: 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?**

# Answer
## 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

```
> 1 votes

# Answer
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.
> 1 votes

---
Tags: r, sequence-alignment, bioconductor
---


# Question
Title: When to account for the blacklisted genomic regions in ChIP-seq data analyses?
We have heard in the group that it is important to keep track of and to filter artifact regions when analysing data from functional genomics experiments, especially ChIP-seq.

Here, we have seen pipelines that remove the ENCODE tracks i) before cross-correlations QCs, ii) after cross-correlations QC but before peak calling and iii) after peak calling.

We have noticed that removal of the tracks does not affect significantly cross-correlation and peak-independent QCs. However, we are not sure whether peak calling should be done on the filtered tracks or not?

# Answer
Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value.

Regardless, if you're using SPP/phantomPeakQual for cross-correlation then note that it already removes the highest peaks from your dataset before computing the cross-correlation (in fact, it can remove most of the actual peaks too, which makes one further wonder what it's actually telling you). I don't know that this is actually documented anywhere, it's something I noticed when going through the code while pondering whether to implement it in deepTools. But at least it's ignoring these regions already :)

In general, it's most convenient to just remove peaks overlapping blacklisted regions. In an ideal world you'd filter out the blacklisted reads before peak calling, but (1) this is really inconvenient (more time and disk required) and (2) I've never seen an appreciable gain in peak calling performance. In theory at least you should be losing sensitivity right around blacklisted regions if you don't remove reads overlapping blacklisted regions, but you have to ask yourself whether you want to trust such peaks anyway. For other QC steps, at least with deepTools we provide a parameter with every tool to specify a BED file of blacklisted regions to skip.

As an aside, there are many fewer blacklisted regions in more recent genome builds (GRCh38 and GRCm38 at least), so this is less of an issue in general with them.
> 10 votes

# Answer
I filtered my bams out of curiosity. I had around 500 peaks before filtering. When i called on the filtered bams, i lost 100 (from the blacklist) and gained 38. The peak call was using macs2. Looking at the gained peaks, they are legit and a binding partner of my TF had peaks there. Note however i ran out of ram for bam \> 8GB using bedtools intersect to filter. And if you do it, sort the bam and make sure the chromosome names are ordered using the same tool between blacklist bed and your bams.
> 0 votes

---
Tags: chip-seq
---


# Question
Title: How to add bootstrap values to the phylogenetic tree generated by OrthoFinder?
When we run the OrthoFinder analysis tool on a group of genomes to get the orthologues shared by them one of the output files include a folder named 'Species\_Tree' that contains a text file named 'SpeciesTree\_rooted.txt' . These are its contents (This is an example I have used for the question - 11 genomes of the *Metabacillus* genus) :

> ((Mlacus:0.244776,Mmangrovi:0.284609)1:0.0277466,(((Mdong:0.0284157,Midriensis:0.0250137)1:0.0553715,Mindicus:0.0904605)1:0.148105,(Mfastidiosus:0.207986,(Mniabensis:0.142138,((Mcrassostreae:0.0124653,Mlitoralis:0.0139468)1:0.116616,(Mhalosaccharovorans:0.0408989,Mschmidt:0.0422379)1:0.0872204)1:0.0402912)1:0.0849362)1:0.0750784)1:0.0277466);

I used iTOL to construct the tree using this file but I could only see the branch lengths and no way to add any bootstrap values if there had been any, to the tree as seen in most research papers - 99, 98, 100 etc. Any suggestions?

# Answer
> 1 votes
Bootstrapping is easy just use **iqtree**

```
iqtree -s file -m GTR+I+G -T 6 -B 1000

```

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.

The `file` is aligned fasta format (its all data pipeline so I forget if I've shifted the format) - although this can be changed to any input format

The result will be in the treefile and is a professional quality tree, complete with bootstraps.

The following up analysis is a Bayesian tree and is more complicated.

**Note** You must resolve/delete gaps using iqtree because of the way the authors are calling gaps - this related to your question here - How to use Gblocks for trimming single-copy gene sequences? . Basically columns with gaps should be deleted.

---

Question, are you sure this is not a Bayesian analysis? It would make sense if it was a Bayesian analysis that is 100% robust (there are robustness scores there). This would make sense for this tree and this data. As far I am aware iTOL is used for viewing and manipulating trees rather than creating them.

---
Tags: phylogenetics, sequence-alignment, phylogeny, taxonomy, orthofinder
---


title: Evolutionary tree from SNP data

text:
# Question
Title: Evolutionary tree from SNP data
I would like to make an evolutionary tree using the SNP data from DArTseq. Are there any software or packages available for this purpose? I came across the software SNphylo, but its home page is not updated. Your suggestions would be appreciated.

# Answer
> 1 votes
At present I would simply use iqtree

```
iqtree -s file -m GTR+I+G -T 6 -B 1000

```

This will parallelise across 6 CPUs and results in 1000 bootstraps, whilst GTR is the mutation model, which models invariant sites under a 4 category gamma distribution. This is a robust model.

You can view the tree in FigTree.

The DNA file is in fasta format

```
> seq1
AGACGT
>seq2
AGAAAT

```

If there's an 'unknown' format bug let me know because I only work in data pipeline so I forget if I've shifted the format.

The result will be in the treefile and is a professional quality tree, complete with bootstraps.

The `I` parameter here could omitted for this dataset because I assume all none of the sites are invariant.

---
Tags: phylogenetics, phylogeny
---


# Question
Title: How to use Gblocks for trimming single-copy gene sequences?
I have run OrthoFinder on a set of 11 genomes and got the results in a folder. From the output folder, I saved all the single-copy gene sequences in a single FASTA file and now wish to remove all gaps using **Gblocks** before using them for constructing a phylogenomic tree. But I don't know the command line used for running **Gblocks**. Any help will be appreciated!

# Answer
Basically this is 16S data. If its proteins there's a different strategy - even if the data set needed is nucleotides (separate question). I would simply use trimal which is available on **conda**

```
trimal -in <inputfile> -out <outputfile> -nogaps

```

The input file is aligned fasta here and will output as aligned fasta.

Any site (alignment position) with an indel `-` will be removed. This is definitely preferred if performing an **iqtree** analysis, because of the way it handles indels. This will also clean up poorly aligned, alignments.

e.g.

```
---ATG
---A-G
--TATG

```

This will result in...

```
AG
AG
AG

```

Ultimately for alignments - there's no magic bullet - and manual editing and inspection is expected e.g. using **seqotron**. There are very good reasons for this concerning circular mathematics. However, modern aligners such as **muscle5** are extremely good.

Trimal is available via conda

```
conda install -c bioconda trimal

```

In terms of using GBlock I'd need to know whether it was simply deleting indel sites, deleting indel sites at \>50% frequency or whether it was performing a genetic diversity calculation per site. I don't think it is doing the latter. It is better to perform diversity saturation calculation explicitly rather than using an unknown automate because thats kinda important. If its deleting sites with excessive \>50% indels - honestly its better to do a blockwise deletion rather than do something half-half. Within SSU there is a reason might be being done because there can be excessive indels at deep taxonomic levels with this data type. Again using iqtree just do it 'blockwise' which is the trimal output.

*Note* I've not read the original GBlock paper, but if you summarise its application I could consider it against **trimal**. However, I definitely prefer **trimal** considering the downstream application.
> 1 votes

---
Tags: phylogenetics, phylogeny, linux, orthologues, orthofinder
---


# Question
Title: Creating a scatterplot by a json data in js
I have a json file of protein data of mouse and human like

```
{
    "name": "P04202",
    "full_name": "Transforming growth factor beta-1 proprotein",
    "symbol": "TGFB1_MOUSE",
    "length": 390,
    "species": "mouse",
    "function": "Transforming growth factor beta-1 proprotein: Precursor of the Latency-associated peptide (LAP) and Transforming growth factor beta-1 (TGF-beta-1) chains, which constitute the regulatory and active subunit of TGF-beta-1, respectively.",
    "nr_modifications": 2821,
    "residues": [
        {
            "pos": 1,
            "residue": "M",
            "modifications": []
        },
        {
            "pos": 2,
            "residue": "P",
            "modifications": []
        },
        ...
        {
            "pos": 60,
            "residue": "A",
            "modifications": [
                {
                    "mod": "[127] Fluoro",
                    "type": "Chemical derivative",
                    "prob": 1
                },
                {
                    "mod": "[53] HNE",
                    "type": "Post-translational",
                    "prob": 0.25
                }
            ]
        },
        ...
        ]
}

```

How do I create a scatter plot \[!\[enter image description here\]\[1\]\]\[1\]

Thanks for any idea

So far I gathered a code to fetch the data from an online json

```
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins()
  }
}

```

This code generates a scatter plot but I will need to hover around the dots

```
script>
    import { extent } from 'd3-array';
    import { scaleLinear } from 'd3-scale';
    export let data = [];

    const margin = 30;

    $: xScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.proteins.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);
</script>

<svg width="800" height="400">
    <g transform="translate({margin}, {margin})">
        {#each data.proteins as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class:human={datapoint.species === 'human'}
                r="20"
            />
        {/each}
        <circle
                cx={10}
                cy={10}
                class="human"
                r="10"
            />
            <circle
            cx={10}
            cy={40}
            class="mouse"
            r="10"
        />

    </g>
</svg>

<style>
    svg {
        border: 1px;
        border-style: solid;
    }
    circle {
        fill: steelblue;
        fill-opacity: 0.5;
    }
    .human {
        fill: red;
    }
</style>

```

# Answer
> 3 votes
It seems like you are already on the right track. The code you provided fetches the protein data from the JSON file and generates a scatter plot with D3.js based on two properties: `length` and `nr_modifications`.

To add an interactive hover event, you can add a tooltip to your scatterplot. The tooltip will display additional information when you hover over each dot.

Below is an example of how you can add a tooltip:

First, add a `<div>` element in your HTML that will be used as the tooltip:

```
<div id="tooltip" style="position: absolute; opacity: 0;">
    <p><strong>Protein:</strong> <span id="protein"></span></p>
    <p><strong>Modifications:</strong> <span id="modifications"></span></p>
    <p><strong>Length:</strong> <span id="length"></span></p>
</div>

```

Then, in your script, add the code to handle the mouseover and mouseout events for each circle:

```
{#each data.proteins as datapoint}
    <circle
        cx={xScale(datapoint.length)}
        cy={yScale(datapoint.nr_modifications)}
        class:human={datapoint.species === 'human'}
        r="20"
        on:mouseover="{(e) => handleMouseover(e, datapoint)}"
        on:mouseout="{(e) => handleMouseout(e)}"
    />
{/each}

```

Define the `handleMouseover` and `handleMouseout` functions in your script. These functions will show and hide the tooltip respectively. For the `handleMouseover` function, it also updates the tooltip content with the protein data:

```
<script>
    // ... your existing code ...

    const handleMouseover = (event, datapoint) => {
        // Update the tooltip content
        document.getElementById('protein').textContent = datapoint.symbol;
        document.getElementById('modifications').textContent = datapoint.nr_modifications;
        document.getElementById('length').textContent = datapoint.length;

        // Show the tooltip
        document.getElementById('tooltip').style.opacity = 1;

        // Position the tooltip
        document.getElementById('tooltip').style.left = event.pageX + 'px';
        document.getElementById('tooltip').style.top = event.pageY + 'px';
    };

    const handleMouseout = () => {
        // Hide the tooltip
        document.getElementById('tooltip').style.opacity = 0;
    };
</script>

```

In standard HTML or JavaScript, the {#each} syntax and usage of {} for expressions and variable binding won't work as they are specific to Svelte. Same applies to the usage of on:mouseover and on:mouseout. In standard HTML or JavaScript, we use onmouseover and onmouseout.

Here are some things which I haven't checked, but you might want to consider looking at:

* It looks like you are trying to load JSON data from an external source, but the function fetchProteins() is never called and its result (the promise returned) is not being awaited. The way it's structured now, Proteins in the returned object will hold a Promise, not the actual data.
* The code assumes that the returned JSON data is an array stored in data.proteins, but the provided JSON example suggests that the data structure might be different. It doesn't show any 'proteins' key in the JSON. The code needs to align with the actual structure of the returned JSON data. `data.proteins` implies that `data` is an object and `proteins` is a property of that object containing an array.
* In your SVG transform attribute, you are using braces {} though the correct syntax for a translate transformation in SVG is translate(x,y).

You might want to try looking at the following changes:

```
export const load = ({ fetch }) => {
  const fetchProteins = async () => {
    const res = await fetch('https:address/assets/ptm.json')
    const data = await res.json()
    return data
  }

  return {
    Proteins: fetchProteins
  }
}

```

```
  $: xScale = scaleLinear()
        .domain(extent(data.map((d) => d.length)))
        .range([0, 800 - (margin * 2)]);

    $: yScale = scaleLinear()
        .domain(extent(data.map((d) => d.nr_modifications)))
        .range([0, 400 - (margin * 2)]);

```

```
<g transform="translate({margin}px, {margin}px)">
        {#each data as datapoint}
            <circle
                cx={xScale(datapoint.length)}
                cy={yScale(datapoint.nr_modifications)}
                class="{datapoint.species === 'human' ? 'human' : ''}"
                r="20"
            />
        {/each}
        <!-- The rest of your code -->
</g>

```

---
Tags: proteins, graphs, json
---


title: How to filter a SAM file by a bed file?

text:
# Question
Title: How to filter a SAM file by a bed file?
I have a large BAM file and I want to filter it by just what's defined in a bed file, and write it to a new file.

Is it possible to do this with samtools, and if so, how? Is there a better tool I could use to do this?

# Answer
Thanks to the helpful comments I figured it out:

```
samtools view -b -h -L bedfile.bed originalbam.bam > newbam.bam

```

-b = Output as bam file
-h = include header
-L = only output alignments overlapping the input BED file
> 5 votes

# Answer
You can also do this with `bedtools intersect`:

```
bedtools intersect -abam input.bam -b bedfile.bed > output.bam

```

This should be identical to the `samtools view` answer. One further feature though is you can output all reads that **don't** overlap with the regions in `bedfile.bed` by adding the `-v` flag.
> 2 votes

---
Tags: sam, samtools, pysam
---


# Question
Title: Getting VCF file that contain common SNPs from 6 VCF file using isec
I have 6 VCF files, where I would like to obtain the SNPs that are common (by position) in all the 6 files. I have tried this command

```
bcftools isec  -p trial2 -n=6 -c all A2_SNP.vcf.gz A_SNP.vcf.gz B_SNP.vcf.gz C_SNP.vcf.gzD_SNP.vcf.gz E_SNP.vcf.gz

```

But this produces 6 files, and through the README file, I think each of the VCF contains the unique SNPs for each of the 6 files, and I could only get sites.txt file that ocntain the positions that are common. Are there any way to get a VCF file that have the common SNPs ?

Thanks

# Answer
> 2 votes
You've used the right tool and gotten the right results. What you need to understand is that VCF contains variant loci as well as the genotypes for \>=1 sample(s) at those loci. The intersect operation will extract common loci, but it does not know which genotypes you're interested in, which is why it generates one file per input, giving you the genotypes from each VCF for the common loci.

You could use the sites file to do a `bcftools merge` operation, limiting to just the common sites using the `-R <regions-file>` option so you'd get a VCF with all the samples across all 6 input VCFs but only for the variant loci.

---
Tags: rna-seq, snp, variant-calling, rna, gatk
---


# Question
Title: Should genotype imputation be ancestry specific?
I'm wondering if imputation, specifically Beagle, needs a reference panel that matches the sample's ancestry group. For example, Beagle documentation suggests the 1000 Genomes Project phase 3 reference panel, which contains all ancestry groups. Is this acceptable to use to impute the genome of, say, an East Asian, African, and European?

# Answer
Basically, yes, you will get better results if you use a reference panel which contains samples with similar ancestry to your targets. Imputation works by matching haplotypes in your target samples to those in a reference panel. If the haplotypes in the target don’t exist in the reference panel, the variants won’t be imputed accurately.

> We demonstrated that using TOPMed sequencing data as the imputation reference panel improves genotype imputation quality in these populations\[Hispanic/Latino and African\], which subsequently enhanced gene-mapping power for complex traits. For rare variants with minor allele frequency (MAF) \< 0.5%, we observed a 2.3- to 6.1-fold increase in the number of well-imputed variants, with 11-34% improvement in average imputation quality, compared to the state-of-the-art 1000 Genomes Project Phase 3 and Haplotype Reference Consortium reference panels.

From this paper.
> 1 votes

# Answer
in my experience you should use a publicly cited reference genome like GRCh37 or GRCh38

https://www.sciencedirect.com/science/article/pii/S0888754317300058

https://gatk.broadinstitute.org/hc/en-us/articles/360035890951-Human-genome-reference-builds-GRCh38-or-hg38-b37-hg19
> 0 votes

---
Tags: imputation, 1000genomes
---


# Question
Title: How can I delete these lines in my fastq file?
Hi guys =) I'm sorry if this is a repeat but I haven't been able to find an answer using the search field or google.

I am trying to edit a fastq file that has a corrupt line. The read in question is VH00362:5:AAAV5KTHV:1:1210:64850:6869 where the read is just a massive repeat – not the expected 110bp. I am trying to remove two lines from this file and asking around I've been told to use `sed` as in

`gzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq.gz`

but when I run `wc -l`  to check the number of lines the results are incoherent, as if it has deleted much more than the two lines I asked. Another issue that I'm seeing with this file is that every time I run `gzip` I get a `invalid compressed data--format violated` error, which tells me the file is corrupted but I'm doing my best to salvage whatever I can from the file as it is an irreplaceable sample.

Is there another way to simply delete these two lines? Any tool or command you could recommend? Any help is deeply appreciated.

# Answer
This is not the thing I do these days but ...

```
gunzip -c file.fastq.gz | sed "45500973,45500974d" > new.fastq

```

You were using `gzip` instead of `gunzip`, and that the source of the error. Finally the out file isn't compressed, so you probably want:

```
gunzip -c file.fastq.gz | sed "45500973,45500974d" | gzip > new.fastq.gz

```

That should work (I tested it and yes it works).
> 2 votes

# Answer
A `fastq` file has four line-separated fields per sequence. If you suspect that a sequence is corrupted, then you would want to delete all four lines related to that sequence. If you just delete the first two lines (`sequence identifier` and `raw sequence letters`), the file will most probably throw an error down the road. In your case, the sequence would look something like this -

```
@VH00362:5:AAAV5KTHV:1:1210:64850:6869
GATTTGGGGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

```

So, the code to delete these four lines from a `fastq.gz` file would be -

```
gunzip -c file.fastq.gz | sed "45500973,45500976d" | gzip > new.fastq.gz

```

This `sed` command specifies a range of lines and says that delete lines from line number `45500973` to `45500976`.
> 1 votes

---
Tags: rna-seq, data-preprocessing
---


# Question
Title: What do the numbers mean in these RNA-Seq gene/transcript TPM files?
From the link https://gtexportal.org/home/datasets, under V7, I'm trying to do R/Python analyses on the Gene TPM and Transcript TPM files. But in these files (and to open them I had to use Universal Viewer since the files are too large to view with an app like NotePad), I'm seeing a bunch of ID's for samples (i.e. `GTEX-1117F-0226-SM-5GZZ7`), followed by transcript ID's like ENSG00000223972.4, and then a bunch of numbers like 0.02865 (and they take up like 99% of the large files). Can someone help me decipher what the numbers mean, please? And are the numbers supposed to be assigned to a specific sample ID? (The amount of letters far exceed the amount of samples, btw). I tried opening these files as tables in R but I do not think R is categorizing the contents of the file correctly.

For context, I am planning to match males with females for sex comparison but in order to do that, I need to get R to categorize everything correctly. (I know that females have "F" where `####-#####-####-#x-#####` where `x` is and males have "M").

**Example subsets**

For gene TPM:

```
#1.2 56202 11688
Name Description GTEX-1117F-0226-SM-5GZZ7 GTEX-111CU-1826-SM-5GZYN GTEX-111FC-0226-SM-5N9B8 GTEX-111VG-2326-SM-5N9BK GTEX-111YS-2426-SM-5GZZQ GTEX-1122O-2026-SM-5NQ91 GTEX-1128S-2126-SM-5H12U GTEX-113IC-0226-SM-5HL5C 7C-0005-SM-57WDL GTEX-ZVE2-0006-SM-51MRW GTEX-ZVZQ-0006-SM-51MR8 GTEX-ZXES-0005-SM-57WCB GTEX-ZXG5-0005-SM-57WCN
ENSG00000223972.4 DDX11L1 0.1082 0.1158 0.02104 0.02329 0 0.04641 0.03076 0.09358 0.121 0.02859 0 0.02373 0.08024 0.02231 0.0189 0.09309 0.07772 0.1596 0.1173 0.07969 0.07763 0.02062 0.07831 0.05159 0.1181 0.04187

```

For transcript TPM:

```
transcript_id gene_id GTEX-1117F-0226-SM-5GZZ7 GTEX-1117F-0426-SM-5EGHI GTEX-1117F-0526-SM-5EGHJ GTEX-1117F-0626-SM-5N9CS GTEX-1117F-0726-SM-5GIEN GTEX-1117F-1326-SM-5EGHH GTEX-1117F-2226-SM-5N9CH GTEX-1117F-2426-SM-5EGGH GTEX-111CU-0226-SM-5GZXC GTEX-111CU-0326-SM-5GZXO
ENST00000373020.4 ENSG00000000003.10 26.84 4.13 13.54 31.17 6.87 39.75 32.25 30.81 19.69 46.42 9.07 27.83 9.04 12.87 18.18 11.12 11.28 7.86 5.52 24.67 6.48 8.13 21.41 5.39 7.98 8.56 51.61 12.51 70.01 17.48 2.1 14.16 3.96 11.39 1.29 6.5 29.93 4.74 3.1 20.82 16.81 14.82 50.8 16.9 9.35 39.3 3.61 2.87 11.73

```

# Answer
The files have pretty self-explanatory headers.

The first file has EnsEMBL gene ID, HGNC symbol and then TPM for each sample. The second file has EnsEMBL gene ID, EnsEMBL transcript ID and TPM for each sample. Each header beginning with `GTEX` is a sample.

You should be using linux command line to look at this, not basic Windows utilities like Notepad. At least use Notepad++, basic Windows utilities are not equipped to handle large text files.
> 1 votes

---
Tags: r, rna-seq, sequence-analysis, genomics, tpm
---


# Question
Title: snakemake restricting resources for specific rule
I am trying to generate several snakemake rules. However, I would like to include some restrictions due to storage limitations

I first would like to copy the the bamfiles, 1 for each sample, to a staging area. However, while I have over 50 samples, I can only process three samples at a time.

Then I want to compare the bam file to chunked candidate files, 1 for each chromosome. This however, I would like to run on 23 cores per sample so that once each bam file is copied 69 cores are working.

finally, if all chunks for a given sample have been processed I want to clean up the staging area, copy the bamfile for the next sample and process the chunks for that sample.

My current implementation, see below, just continues copying bamfiles and crashes when I run out of strorage;

```
rule copy_bam:
    input:
        bamfile=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        bamindex=os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai")
    output:
        staging_bamfile="{staging_area}/{sample}.bam",
        staging_bamindex="{staging_area}/{sample}.bam.bai"
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bamfile} {output.staging_bamfile}
        cp {input.bamindex} {output.staging_bamindex}
        """

rule process_chunks:
    input:
        chunkfile="{staging_area}/chunks/{chunk}_with_header.bcf.gz",
        bamfile=rules.copy_bam.output.staging_bamfile
    output:
        dynamic("{staging_area}/calls/{sample}/{chunk}.bcf")

    threads: 23
    shell:
        """
        process chunks
        """

rule cleanup:
    input:
        bamfile="{staging_area}/{sample}.bam",
        bamindex="{staging_area}/{sample}.bam.bai",
        processed=lambda wildcards: expand("{staging_area}/calls/{sample}/{chunk}.bcf", chunk=get_chunks()[wildcards.sample], staging_area=config['staging_area'], sample=wildcards.sample)
    output:
        touch("{staging_area}/{sample}/cleanup_done.txt")
    shell:
        """
        rm {input.bamfile}
        rm {input.bamindex}
        rm {input.processed}
        """

```

run with: `snakemake --resources bam_copy=3 --cores 70`

# Answer
One Snakemake feature you can use is temporary outputs.

If you make your outputs from the `copy_bam` rule temporary, they'll be cleared up after no other rule requires that file as an input. Your Snakefile would look something like this:

```
rule copy_bam:
    input:
        bam = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam"),
        idx = os.path.join(config['bam_folder'], "{sample}.gatk.sorted.MD.bam.bai"),
    output:
        bam = temp("{staging_area}/{sample}.bam"),
        idx = temp("{staging_area}/{sample}.bam.bai"),
    wildcard_constraints:
        sample="[A-Za-z0-9_]+"
    resources:
        bam_copy=1
    shell:
        """
        cp {input.bam} {output.bam}
        cp {input.idx} {output.idx}
        """
# omitted for brevity

```

You can then specify how many jobs run concurrently with the `--cores`/`--jobs` options to keep the space usage small.
> 2 votes

---
Tags: snakemake
---


# Question
Title: StringTie discovers less novel isoforms when reference annotation is provided
I am analyzing PacBio IsoSeq data and I am using StringTie to assemble transcripts. I noticed the difference in the output files when I add a reference genome annotation (-G option in StringTie). In general I have less transcripts assembled when I give as imput reference genome annotation and I was wondering why. I know that the annotation of the genome is not the best (Borelian orangutang) and I would rather trust the data rather than the annotation provided.

I am attaching an example for one of the genes (IGV): the first track is reference annotation; the second is StringTie output when no reference annotation is provided; the third is StringTie output when reference annotation is provided (in the third case two novel assembled isoforms are missing)

# Answer
Not very familiar with what StringTie is doing, but in general reference annotations are helpful in providing a prior for what the tool expects.

I would suggest looking at the manual and using some of the provided options and outputs to probe the coverage of the various transcript species (I believe these are in the "attributes" section of the GTF output).

If the transcripts that are missing when you use `-G` have a coverage of just a couple of reads, they are likely to be less interesting. On the other hand, if they have relatively high coverage, then you have more reason to be skeptical of the existing annotation.

My guess is that the tool is more heavily weighting existing transcripts based on the annotation, and therefore excluding the extra transcripts as being below some threshold.

I'd recommend consulting the "recommended workflow" in the manual that supposedly allows for novel isoforms. Ensure that you are using this workflow and not the more streamlined workflow, which does not allow for novel isoforms.

See also the comment suggesting comparing to other IsoSeq-oriented tools.
> 1 votes

---
Tags: pacbio, isoform, stringtie
---


# Question
Title: Easier / lazy way to convert an existing bash pipeline to nextflow pipeline?
Not going to lie, my head hurts when I read the NextFlow documentation (No dig on the documentation, I think my head is losing its processing capacity lol).

I had written a pipeline in Bash, essentially iterates over a text file containing full paths of R1 and R2 reads, followed by a series of commands and liberal use of gnu parallel to improve efficiency.

Is there any template I can follow, where I can slot in my commands and get things running quickly?

EDIT - The skeleton of pipeline goes as follows (As I'm not sure if I can put our inhouse pipeline in its entirety here) -

Bash code which loops through R1 and R2, aligns with BWA-mem2 and sorts the bam, then using gnu parallel, I essentially follow the steps as outlined in GATK Best Practices (looking to drop BQSR, but the steps remain) to obtain a gvcf file. The steps in Best Practices have been turned to functions, rather than each step being called sequentially.

# Answer
> 5 votes
If your objective is to convert Bash scripts into Nextflow code, you might find Viash (https://viash.io) a particularly useful open source tool. It automatically generates independent Nextflow modules from each of your scripts, along with the corresponding Docker containers for each module. Docker containers are generated on the fly based on your config description in Viash. These modules can be coupled together to build a (modular) Nextflow pipeline.

# Answer
> 2 votes
The gist would be something like:

```
fastqs = Channel.of(...some annoying groovy code to make this a list of tuples...)
workflow{
  align(fastqs)
  makeGVCF(align.out.bams)
  ... collect() the output from above joint call them all...
}

process align {
  input:
    tuple val(sampleName), path(reads)
  output:
    tuple val(sampleName), path("${sampleName}.bam"), emit: bams
  """
  bwa ... | samtools sort ...
  """
}

process makeGVCF{
  input:
    tuple val(sampleName), path(bam)
  output:
    path "${sampleName}.gvcf", emit: gvcfs
  """
  ...haplotype caller code...
  """
}

```

# Answer
> 2 votes
Since you asked about templates, if your task is something rather generic, there is a good chance that a module, or several, that you can use to build your workflow exists in https://nf-co.re/modules Generally I would recommend taking a look at the source of nf-core pipelines (e.g. https://nf-co.re/sarek for variant calling) to get an idea how to get your samplesheets in.

---
Tags: bash, nextflow
---


title: Singularity And CNV Kit

text:
# Question
Title: Singularity And CNV Kit
I would like to try an use Singularity with CNV kit as the HPC I use seems to work better with Singularity:

***Initial Attempt:***

```
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "ls /"

```
* *Result*:

```
bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var

```

***Attempt 1:***

```
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /proc/ && ls"

```

This results in a folder with a series of numbers but not the cnvkit py files.

***Attempt 2:***

```
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /singularity/ && ls"

```

This results in `/bin/bash: line 1: cd: /singularity/: Not a directory`

***Attempt 3:***

```
 singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "cd /media/ && ls"

```

This results in nothing being printed to the screen. I tried this command three times to be sure.

How can I find the Singularity path I would need to use with `singularity exec` to run from the Singularity image.

# Answer
I like this question a lot. Firstly I don't know singularity, but I know Docker and it will be similar. `docker exec` is very similar.

So for example the reasons are likely as follows:

* `media` could be empty and hence no output.
* the Python script may not be in `/proc`
* singularity is likely to be a file rather than a directory.

One solution is,

```
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "ls -l /"

```

i.e. ask it which are the directories. Another solution is..

```
singularity exec ~/CNV_kit/cnvkit_latest.sif /bin/bash -c "find / -name 'cnvkit*'"

```

This will locate `cnvkit.py`, note the use of ' and " for this sort of command.

Personally, I would dive into the container and issue the command of interest, manually. This way you get a feel for how it's constructed and don't need to keep issuing commands from outside the container. In Docker this is easy and it's impossible to keep a third party out of a container (even if the developer tries). In singularity it is likely similar.

---

**From the notes** I would simply access the singularity container directly.

Accessing a singularity container is different from Docker (I looked it up) ...

```
singularity shell ~/CNV_kit/cnvkit_latest.sif

```

Once inside you'll be able to look around, try out the `find` command, look in the `bin` etc... There will be no need to issue commands from outside the container to get it working.

---

From further comments,

A node on an HPC is a group of usually 8 to 32+ CPUs. What you have is a 'master' or 'head' node, where you log into and keep your data, and the stuff you do the work on is called 'nodes'. It has a more specific name, but it's not very cool, so we don't use it. The HPC will comprise a large number of nodes, e.g. 100+, but can be smaller particularly with the massive CPU machines around now. A node concept is a bit different on a NVidia machine if its part of a cluster.

The issue is someone could have called a script 'node', which would be confusing. It's complicated, because someone could have their own research machine which is attached and maintained on. the HPC, but only they can access and called it e.g. 'Joes\_node', of 'Jills\_node'.

Personally, I would sit down with your administrator and ask them to explain singularity and this the comments herein and this. They will know that at admin level (thats their job).

From my point of view, if I'd access to a cluster using singularity I could figure it out. I use GCP now, it's geeky but cool, so my information is out of date on non-Cloud clusters.

---

**RamRS comments** What RamRS is saying and what I'm trying to assess whether you are within a singularity container. Again, I don't know the architecture.

**Merging With Accepted Answer:**

The main file in cnvkit is `cnvkit.py` in order to find that you need to look using `singularity exec`. The following command gives the path to the cnvkit.py file:

```
singularity exec cnvkit_latest.sif /bin/bash -c "cd /usr/local/bin/ && ls"

```

An example of how to use this path is below:

```
# Directory containing sample files:
directory="/path/to/bam/files/folder"

# Define path to singularity image:
cnvkit_container="/path/to/cnvkit_latest.sif"

cnvkit_path="/usr/local/bin/cnvkit.py"

singularity exec "$cnvkit_container" "$cnvkit_path" call "$directory/sample.cns" -o "$directory/sample.call.cns"

```
> 2 votes

# Answer
Thanks to OP for sharing that their image is from https://hub.docker.com/r/etal/cnvkit

I downloaded this and it seems pretty simple.

My tests:

```
singularity pull docker://etal/cnvkit:latest
$ singularity shell ./cnvkit_latest.sif
INFO:    Converting SIF file to temporary sandbox...
bash: warning: setlocale: LC_ALL: cannot change locale (en_US.UTF-8)
Singularity> ls /
bin  boot  dev  environment  etc  home  lib  lib32  lib64  libx32  media  mnt  opt  proc  root  run  sbin  singularity  srv  sys  tmp  usr  var
Singularity> ls /singularity
/singularity
Singularity> cnvkit.py
usage: cnvkit.py [-h]
                 {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                 ...

CNVkit, a command-line toolkit for copy number analysis.

positional arguments:
  {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                        Sub-commands (use with -h for more info)
    batch               Run the complete CNVkit pipeline on one or more BAM files.
    target              Transform bait intervals into targets more suitable for CNVkit.
    access              List the locations of accessible sequence regions in a FASTA file.
    antitarget          Derive off-target ("antitarget") bins from target regions.
    autobin             Quickly calculate reasonable bin sizes from BAM read counts.
    coverage            Calculate coverage in the given regions from BAM read depths.
    reference           Compile a coverage reference from the given files (normal samples).
    fix                 Combine target and antitarget coverages and correct for biases. Adjust raw coverage data according to the given reference, correct potential biases and re-center.
    segment             Infer copy number segments from the given coverage table.
    call                Call copy number variants from segmented log2 ratios.
    diagram             Draw copy number (log2 coverages, segments) on chromosomes as a diagram. If both the raw probes and segments are given, show them side-by-side on each chromosome (segments on the left side, probes on the right
                        side).
    scatter             Plot probe log2 coverages and segmentation calls together.
    heatmap             Plot copy number for multiple samples as a heatmap.
    breaks              List the targeted genes in which a copy number breakpoint occurs.
    genemetrics         Identify targeted genes with copy number gain or loss.
    sex                 Guess samples' sex from the relative coverage of chromosomes X and Y.
    metrics             Compute coverage deviations and other metrics for self-evaluation.
    segmetrics          Compute segment-level metrics from bin-level log2 ratios.
    bintest             Test for single-bin copy number alterations.
    import-picard       Convert Picard CalculateHsMetrics tabular output to CNVkit .cnn files. The input file is generated by the PER_TARGET_COVERAGE option in the CalculateHsMetrics script in Picard tools. If 'antitarget' is in the
                        input filename, the generated output filename will have the suffix '.antitargetcoverage.cnn', otherwise '.targetcoverage.cnn'.
    import-seg          Convert a SEG file to CNVkit .cns files.
    import-theta        Convert THetA output to a BED-like, CNVkit-like tabular format. Equivalently, use the THetA results file to convert CNVkit .cns segments to integer copy number calls.
    import-rna          Convert a cohort of per-gene log2 ratios to CNVkit .cnr format.
    export              Convert CNVkit output files to another format.
    version             Display this program's version.

options:
  -h, --help            show this help message and exit

See the online manual for details: https://cnvkit.readthedocs.io
Singularity> which cnvkit.py
/usr/local/bin/cnvkit.py

```

As you can see, since `cnvkit.py` is in `/usr/local/bin/`, simply typing `cnvkit.py` within the shell works. Let's try an exec with no frills:

```
$ singularity exec cnvkit_latest.sif cnvkit.py
INFO:    Converting SIF file to temporary sandbox...
usage: cnvkit.py [-h]
                 {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                 ...

CNVkit, a command-line toolkit for copy number analysis.

positional arguments:
  {batch,target,access,antitarget,autobin,coverage,reference,fix,segment,call,diagram,scatter,heatmap,breaks,genemetrics,gainloss,sex,gender,metrics,segmetrics,bintest,import-picard,import-seg,import-theta,import-rna,export,version}
                        Sub-commands (use with -h for more info)
    batch               Run the complete CNVkit pipeline on one or more BAM files.
    target              Transform bait intervals into targets more suitable for CNVkit.
    access              List the locations of accessible sequence regions in a FASTA file.
    antitarget          Derive off-target ("antitarget") bins from target regions.
    autobin             Quickly calculate reasonable bin sizes from BAM read counts.
    coverage            Calculate coverage in the given regions from BAM read depths.
    reference           Compile a coverage reference from the given files (normal samples).
    fix                 Combine target and antitarget coverages and correct for biases. Adjust raw coverage data according to the given reference, correct potential biases and re-center.
    segment             Infer copy number segments from the given coverage table.
    call                Call copy number variants from segmented log2 ratios.
    diagram             Draw copy number (log2 coverages, segments) on chromosomes as a diagram. If both the raw probes and segments are given, show them side-by-side on each chromosome (segments on the left side, probes on the right
                        side).
    scatter             Plot probe log2 coverages and segmentation calls together.
    heatmap             Plot copy number for multiple samples as a heatmap.
    breaks              List the targeted genes in which a copy number breakpoint occurs.
    genemetrics         Identify targeted genes with copy number gain or loss.
    sex                 Guess samples' sex from the relative coverage of chromosomes X and Y.
    metrics             Compute coverage deviations and other metrics for self-evaluation.
    segmetrics          Compute segment-level metrics from bin-level log2 ratios.
    bintest             Test for single-bin copy number alterations.
    import-picard       Convert Picard CalculateHsMetrics tabular output to CNVkit .cnn files. The input file is generated by the PER_TARGET_COVERAGE option in the CalculateHsMetrics script in Picard tools. If 'antitarget' is in the
                        input filename, the generated output filename will have the suffix '.antitargetcoverage.cnn', otherwise '.targetcoverage.cnn'.
    import-seg          Convert a SEG file to CNVkit .cns files.
    import-theta        Convert THetA output to a BED-like, CNVkit-like tabular format. Equivalently, use the THetA results file to convert CNVkit .cns segments to integer copy number calls.
    import-rna          Convert a cohort of per-gene log2 ratios to CNVkit .cnr format.
    export              Convert CNVkit output files to another format.
    version             Display this program's version.

options:
  -h, --help            show this help message and exit

See the online manual for details: https://cnvkit.readthedocs.io
INFO:    Cleaning up image...

```

The answer seems to be a lot easier than having to find `cnvkit.py`, as it is already in the PATH.
> 0 votes

---
Tags: bash, singularity
---


# Question
Title: Finding transposable elements using RepeatMasker
I'm using RepeatMasker to detect, classify the Transposable elements. My Input is a eukaryotic non-reference genome. I made a run via RepeatMasker many times to Mask the TEs, but return 0 Annotation tables. Further, I used a different `-species` option each time, but Annotation results contain only short and simple repeats.

Do I have to create a unique library to detect them?

The command I used:

```
./RepeatMasker -no_is -noint - species mammal Genome.fna

```

Results:

```
sequences:         32573
total length: 2004063690 bp  (1981588036 bp excl N/X-runs)
GC level:         41.28 %
bases masked:   33282890 bp ( 1.66 %)
==================================================
               number of      length   percentage
               elements*    occupied  of sequence
--------------------------------------------------
SINEs:                0            0 bp    0.00 %
      ALUs            0            0 bp    0.00 %
      MIRs            0            0 bp    0.00 %

LINEs:                0            0 bp    0.00 %
      LINE1           0            0 bp    0.00 %
      LINE2           0            0 bp    0.00 %
      L3/CR1          0            0 bp    0.00 %

LTR elements:         0            0 bp    0.00 %
      ERVL            0            0 bp    0.00 %
      ERVL-MaLRs      0            0 bp    0.00 %
      ERV_classI      0            0 bp    0.00 %
      ERV_classII     0            0 bp    0.00 %

DNA elements:         0            0 bp    0.00 %
     hAT-Charlie      0            0 bp    0.00 %
     TcMar-Tigger     0            0 bp    0.00 %

Unclassified:         0            0 bp    0.00 %

Total interspersed repeats:        0 bp    0.00 %

Small RNA:            0            0 bp    0.00 %

Satellites:           1          267 bp    0.00 %
Simple repeats:  663679     27953718 bp    1.39 %
Low complexity:  123705      6154714 bp    0.31 %
==================================================

```

# Answer
That's one answer ^, not sure if you NEED to do that, which wouldn't actually classify your elements as LTRs, it would just find repetitive sequences in your new genome (as you discovered).

To get them classified one standard approach is to use evidence of homology to sequences of LTR elements with classifications already, and you can get those in Repbase. If you want to download Repbase, you have to sign up for an account with GIRI but I think they give to anyone who asks. Then you can download the section of Repbase you're interested in (e.g. the whole thing, LTR, DNA-TE, ...) and provide the path to the downloaded library (in FASTA format I think) using `-lib`, like

`RepeatMasker -lib Repbase_LTR.fa genome.fa`

Then your RepeatMasker output table should have info about other kinds of TEs.

If you want to check out my approach to classification by homology you might find some ideas you'd like to try. I used Dfam+nhmmer and Repbase+tblastx to identify evidence for LTR retrotransposon classification.

My software pipeline is available as open-source software. If you want detailed high-quality annotations of LTR retrotransponsons, PhyLTR will automatically annotate putative protein-coding sequence regions, classify them, and remove false-positives like a tandem array of DNA TEs which look like LTR-Rs (good diagram of LTR-R false positives here: https://github.com/oushujun/LTR\_retriever/blob/master/Manual.pdf)

Phylogenetic Analysis of LTR retrotransposons https://github.com/mcsimenc/PhyLTR
> 2 votes

# Answer
You need to model the repeats in your de novo genome.

*See: http://www.repeatmasker.org/RepeatModeler/ and https://www.biostars.org/p/154290/ (from previous answer)*

You first build a database named "name\_of\_your\_database" (insert whatever you want to call it) and run it on the fasta file of your genome (e.g. your\_genome.fasta).

Make you find the right paths for the software.

```
## build the database
BuildDatabase -engine ncbi -name name_of_your_database your_genome.fasta

## run the modeler (in this case using ncbi)
RepeatModeler -database name_of_your_database -engine ncbi

```
> 3 votes

# Answer
In case someone finds this years later: OP's problem was the use of the "-noint" flag. It stands for "no interspersed", so the OP explicitly asked RepeatMasker NOT to look for transposable elements. If they omitted the flag, RepeatMasker would have automatically downloaded mammal TE sequences (nowadays it uses Dfam). Even if they had a custom library like suggested above, they would still need to delete "-noint". Check your commands!
> 1 votes

---
Tags: repeat-elements, transposable-elements
---


# Question
Title: Why does this nextflow script finish after running one sample?
I have 36 samples in total in the bam\_files folder, with the name like this 20230306\_CH\_EP\_C01.md.bam

I expect the code to output all the 36 samples, one by one. But the run stopped after only running one sample. I have searched online, but did not find a solution. I am using nextflow version 2.

```
bam_chl = Channel
     .fromPath("bam_files/*.md.bam")
     .map{ item -> tuple(item.baseName.substring(0, 18), item) }

ref_chl = Channel
      .fromPath("sequences/hg38.fa")
ref_index_chl = Channel
      .fromPath("sequences/hg38.fa.fai")
amplicon_chl = Channel
      .fromPath("sequences/tp53_amplicon.hg38.interval_list")
target_chl = Channel
      .fromPath("sequences/TP53.hg38.interval_list")

process get_pcr_metric {

  publishDir path: 'pcr_metric', pattern: '*.csv'

  input:
  tuple val(pair_id), path(bam)
  path(ref)
  path(ref_index)
  path(amplicon)
  path(target)

  output:
  path("*.pcr_metrics.csv")

  """
  java -jar /picard.jar CollectTargetedPcrMetrics \
       I=$bam \
   O=${pair_id}_pcr_metrics.txt \
       R=$ref \
   AMPLICON_INTERVALS=$amplicon \
       TARGET_INTERVALS=$target

  sed -n -e 7p -e 8p ${pair_id}_pcr_metrics.txt |  sed 's/\t/,/g'   > ${pair_id}.pcr_metrics.csv

  """

}

workflow {

get_pcr_metric(bam_chl, ref_chl,ref_index_chl, amplicon_chl, target_chl)
}

```

# Answer
The issue here is that the `fromPath` factory method produces a queue channel, but you don't want all of your input channels to be queue channels. Most of the time, what you want is one queue channel and one or more value channels when you require multiple input channels.

You *could* use the `value` factory method to create the value channels, but a simpler way is to just pass in some file objects. The latter works because a value channel is implicitly created by a process when it is invoked with a simple value. For example:

```
params.ref_fasta = "./sequences/hg38.fa"
params.amplicon_intervals = "./sequences/tp53_amplicon_intervals.hg38.interval_list"
params.target_intervals = "./sequences/TP53.hg38.interval_list"

params.bam_files = "./bam_files/*.md.bam"
params.outdir = './results'

```

```
process get_pcr_metric {

    publishDir path: "${params.outdir}/pcr_metric", mode: 'copy'

    input:
    tuple val(pair_id), path(bam_file)
    path ref_fasta
    path amplicon_intervals
    path target_intervals

    output:
    path "${pair_id}_pcr_metrics.txt"

    """
    picard CollectTargetedPcrMetrics \\
         I="${bam_file}" \\
         O="${pair_id}_pcr_metrics.txt" \\
         R="${ref_fasta}" \\
         AMPLICON_INTERVALS="${amplicon_intervals}" \\
         TARGET_INTERVALS="${target_intervals}"
    """
}

```

```
workflow {

    ref_fasta = file( params.ref_fasta )
    amplicon_intervals = file( params.amplicon_intervals )
    target_intervals = file( params.target_intervals )

    Channel
        .fromPath( params.bam_files )
        .map { tuple( it.getBaseName(2), it ) }
        .set { bam_files }

    get_pcr_metric( bam_files, ref_fasta, amplicon_intervals, target_intervals )
}

```
> 3 votes

# Answer
The code can be made to work on *all* samples by changing the interface of the process and using the `combine` operator:

```
ref_chl = Channel
      .fromPath("sequences/hg38.fa")
ref_index_chl = Channel
      .fromPath("sequences/hg38.fa.fai")
amplicon_chl = Channel
      .fromPath("sequences/tp53_amplicon.hg38.interval_list")
target_chl = Channel
      .fromPath("sequences/TP53.hg38.interval_list")

bam_chl = Channel
     .fromPath("bam_files/*.md.bam")
     .map{ item -> tuple(item.baseName.substring(0, 18), item) }
     .combine(ref_chl)
     .combine(ref_index_chl)
     .combine(amplicon_chl)
     .combine(target_chl)

process get_pcr_metric {

  publishDir path: 'pcr_metric', pattern: '*.csv'

  input:
  tuple val(pair_id), path(bam), path(ref), path(ref_index), path(amplicon), path(target)

  output:
  path("*.pcr_metrics.csv")

  """
  java -jar /picard.jar CollectTargetedPcrMetrics \
       I=$bam \
   O=${pair_id}_pcr_metrics.txt \
       R=$ref \
   AMPLICON_INTERVALS=$amplicon \
       TARGET_INTERVALS=$target

  sed -n -e 7p -e 8p ${pair_id}_pcr_metrics.txt |  sed 's/\t/,/g'   > ${pair_id}.pcr_metrics.csv
  """
}

workflow {
  get_pcr_metric(bam_chl)
}

```

The `combine` operator as the name implies generates all combinations. The different file objects are then in one tuple.

Since only one `Channel` is passed to the process, it's now also possible to rewrite the workflow code as follows:

```
workflow {

  refs = [
    "sequences/hg38.fa",
    "sequences/hg38.fa.fai",
    "sequences/tp53_amplicon.hg38.interval_list",
    "sequences/TP53.hg38.interval_list"
  ]

  ref_chl = Channel.from(refs)
   | map{ pointer -> file(pointer) }
   | collect

  Channel.fromPath("bam_files/*.md.bam")
   | map{ item -> tuple(item.baseName.substring(0, 18), item) }
   | combine( ref_chl )
   | get_pcr_metric

}

```
> 2 votes

---
Tags: nextflow
---


title: ProteinPromt installation

text:
# Question
Title: ProteinPromt installation
i dont know how to build/pull this container (https://codebase.helmholtz.cloud/proteinprompt/ProteinPrompt) can you advice me how to install it?

EDIT:

My soft - Ubuntu 20.04 lts

Thanks to @M\__ I manage to build image (almost). But I have an issue :

```
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker build -t dockerfile .
[+] Building 1.0s (15/15) FINISHED
 => [internal] load build definition from Dockerfile                       0.0s
 => => transferring dockerfile: 1.19kB                                     0.0s
 => [internal] load .dockerignore                                          0.0s
 => => transferring context: 2B                                            0.0s
 => [internal] load metadata for docker.io/continuumio/miniconda3:4.10.3-  0.9s
 => [ 1/10] FROM docker.io/continuumio/miniconda3:4.10.3-alpine@sha256:e0  0.0s
 => [internal] load build context                                          0.0s
 => => transferring context: 650B                                          0.0s
 => CACHED [ 2/10] RUN rm /bin/sh && ln -s /bin/bash /bin/sh               0.0s
 => CACHED [ 3/10] RUN conda update -n base conda                          0.0s
 => CACHED [ 4/10] COPY ./data/proteinPrompt.yml /                         0.0s
 => CACHED [ 5/10] RUN conda env update --name root -f /proteinPrompt.yml  0.0s
 => CACHED [ 6/10] RUN mkdir /proteinPrompt/                               0.0s
 => CACHED [ 7/10] COPY ./proteinPrompt/*py /proteinPrompt/                0.0s
 => CACHED [ 8/10] RUN mkdir /model/                                       0.0s
 => CACHED [ 9/10] RUN wget https://zenodo.org/record/5708224/files/rf_hu  0.0s
 => CACHED [10/10] COPY ./data/all.proteins.HUMAN.autocorrelation.pbz2 /m  0.0s
 => exporting to image                                                     0.0s
 => => exporting layers                                                    0.0s
 => => writing image sha256:ae93ff408fb6468280a5f0344bb2ac5ac811df2e5d32a  0.0s
 => => naming to docker.io/library/dockerfile                              0.0s
WARNING: buildx: failed to read current commit information with git rev-parse --is-inside-work-tree

```

and image is named 'dockerfile'? and it do not work

```
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker images
REPOSITORY    TAG           IMAGE ID       CREATED          SIZE
dockerfile    latest        ae93ff408fb6   11 minutes ago   2.02GB
alphafold     latest        c41a4b5e9d48   13 days ago      10GB
nvidia/cuda   11.0.3-base   4f3e95122fab   23 months ago    122MB

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker run  -v $(pwd)/data:/data proteinprompt search -f data/input.fa -o data/output.csv
Unable to find image 'proteinprompt:latest' locally
docker: Error response from daemon: pull access denied for proteinprompt, repository does not exist or may require 'docker login': denied: requested access to the resource is denied.
See 'docker run --help'.

```

EDIT 2:

```
I try to edit dockerfile (in several places)

FROM continuumio/miniconda3:4.10.3-alpine
MAINTAINER Sebastian Canzler <sebastian.canzler@ufz.de>
LABEL authors="sebastian.canzler@ufz.de" \
    description="Docker image of ProteinPrompt - a tool to predict protein-protein interactions based on amino acid sequences."
RUN apt-get update && apt-get install -y git #<-----HERE

#Try This to get a bash in order to run 'source activate modfinder'
RUN rm /bin/sh && ln -s /bin/bash /bin/sh

RUN conda update -n base conda
COPY ./data/proteinPrompt.yml /
RUN conda env update --name root -f /proteinPrompt.yml && \
    rm -rf /opt/conda/pkgs/*
ENV PATH /opt/conda/envs/proteinPrompt/bin:$PATH
# copy the proteinPrompt files
RUN mkdir /proteinPrompt/
COPY ./proteinPrompt/*py /proteinPrompt/
ENV PATH /proteinPrompt:$PATH

# get the random forest model from zenodo
RUN mkdir /model/
#COPY ./data/rf_human_750_inv.rev.pkl /model/random_forest.pkl
RUN wget https://zenodo.org/record/5708224/files/rf_human_750_inv.rev.pkl?download=1 -O /model/random_forest.pkl
COPY ./data/all.proteins.HUMAN.autocorrelation.pbz2 /model/all.proteins.HUMAN.autocorrelation.pbz2

# make proteinPrompt executable
ENTRYPOINT ["python", "/proteinPrompt/ProteinPrompt.py"]
#CMD ["-h"]

```

but there is a bash problem:

```
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker build -t dockerfile .
[+] Building 1.6s (10/15)
 => [internal] load .dockerignore                                          0.0s
 => => transferring context: 2B                                            0.0s
 => [internal] load build definition from Dockerfile                       0.0s
 => => transferring dockerfile: 1.23kB                                     0.0s
 => [internal] load metadata for docker.io/continuumio/miniconda3:4.10.3-  1.0s
 => [ 1/11] FROM docker.io/continuumio/miniconda3:4.10.3-alpine@sha256:e0  0.0s
 => [internal] load build context                                          0.0s
 => => transferring context: 650B                                          0.0s
 => CACHED [ 2/11] RUN rm /bin/sh && ln -s /bin/bash /bin/sh               0.0s
 => CACHED [ 3/11] RUN conda update -n base conda                          0.0s
 => CACHED [ 4/11] COPY ./data/proteinPrompt.yml /                         0.0s
 => CACHED [ 5/11] RUN conda env update --name root -f /proteinPrompt.yml  0.0s
 => ERROR [ 6/11] RUN apt-get update && apt-get install -y git             0.5s
------
 > [ 6/11] RUN apt-get update && apt-get install -y git:
#0 0.537 /bin/ash: apt-get: not found
------
WARNING: buildx: failed to read current commit information with git rev-parse --is-inside-work-tree
Dockerfile:15
--------------------
  13 |         rm -rf /opt/conda/pkgs/*
  14 |     ENV PATH /opt/conda/envs/proteinPrompt/bin:$PATH
  15 | >>> RUN apt-get update && apt-get install -y git
  16 |
  17 |     # copy the proteinPrompt files
--------------------
ERROR: failed to solve: process "/bin/ash -eo pipefail -c apt-get update && apt-get install -y git" did not complete successfully: exit code: 127

```

EDIT 3:

```
roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker pull registry.hzdr.de/proteinprompt/proteinprompt:latest
latest: Pulling from proteinprompt/proteinprompt
5843afab3874: Already exists
72c585c386a3: Already exists
87d798ddb41c: Already exists
316f048b0253: Pull complete
ae11eb6b6512: Pull complete
d3dc3767a31c: Pull complete
ea3dbaf26ce1: Pull complete
061c354f8f77: Pull complete
60eb5484e32a: Pull complete
0625cf0fe317: Pull complete
442818fdc663: Pull complete
203550cafc1f: Pull complete
Digest: sha256:01f0251f34313437e5fcc2dbb26d1ec1103d50aef6335b2261a2048be8ed0813
Status: Downloaded newer image for registry.hzdr.de/proteinprompt/proteinprompt:latest
registry.hzdr.de/proteinprompt/proteinprompt:latest

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker images
REPOSITORY                                     TAG           IMAGE ID       CREATED         SIZE
alphafold                                      latest        c41a4b5e9d48   13 days ago     10GB
registry.hzdr.de/proteinprompt/proteinprompt   latest        cd99d22f3cb2   18 months ago   1.93GB
nvidia/cuda                                    11.0.3-base   4f3e95122fab   23 months ago   122MB

roboczy@marcin-B550-AORUS-ELITE-AX-V2:/home/marcin/ProteinPrompt$ docker run  -v $(pwd)/data:/data proteinprompt search -f data/input.fa -o data/output.csv
Unable to find image 'proteinprompt:latest' locally
docker: Error response from daemon: pull access denied for proteinprompt, repository does not exist or may require 'docker login': denied: requested access to the resource is denied.
See 'docker run --help'.

```

Do you know what is wrong?

---

Final approach,

```
docker pull registry.hzdr.de/proteinprompt/proteinprompt:latest

```

# Answer
The docker file is here https://codebase.helmholtz.cloud/proteinprompt/ProteinPrompt/-/blob/main/Dockerfile

You need to download and install a copy of Docker (below).

These are the instructions how to build the container locally from the remote Docker file. Copy and paste this code into a file called `Dockerfile`. Here I placed it in a `tmp` directory, then simply (see **note**)

```
 docker build -t tmp/dockerfile .

```

Just to note the authors may have specific instructions for the build command.

**How to install Docker** Before this command is run you need to install docker here https://www.docker.com/products/docker-desktop/

After installation you need a special command to run Docker in sudo-less mode is:

```
sudo usermod -aG docker

```

The first command I gave will then build the Dockerfile. There quite a lot of commands BTW to run/execute a Dockerfile, but you can also ask more questions.

---

*Why am I recommending this way*?

Simple, because normally downloading docker is done via https://docker.com and you would issue a `pull` request. For some weird reason the authors haven't done this and I don't know what the outcome of an external download is. If it is built via the Dockerfile then Docker will recognise the Docker container. Once the container is built locally you can upload it to your own Docker site and delete it locally and access it any time via a `pull` request or a `run` request (will automate a `pull` request).

Docker is very cool BTW, worth learning.

---

**Comments**, following bug report

I don't think you've got git installed (`apt-get install -y git`). Once that is installed I would would try again. **BUT FIRST** I would delete the current build docker rmi dockerfile. If docker rmi ... fails (I don't think it will) then remove all containers before deleting the image with `docker rm $(docker ps -qa)`. The use `rmi` and then install git `apt-get install -y git` and do the docker build command.

If this doesn't work (same error) ..... you need to alter the Dockerfile and build it again from the new file .. Its a bug by the authors simply add the line `RUN apt-get update && apt-get install -y git` ... This is because there's a chance that git needs to be on the inside of the container not the outside.

If all the above fails, I would download the entire image/container from the website the authors have loaded it. I don't know how to pull from a non-docker website, I've never done it. Note I'll place these comments in an updated post. I am honestly skeptical about the authors ability to write/develop Docker.
> 1 votes

---
Tags: software-installation, protein-protein-interaction
---


title: Integrate bulk RNA and ATAC-seq genes

text:
# Question
Title: Integrate bulk RNA and ATAC-seq genes
I have up and down-regulated genes from bulk RNA result from DeSeq2. Genes that are differential in chromatin accessibility between two conditions from DiffBind. I would like to find the genes overlap between bulk RNA and ATAC-seq assay as these papers but don't know the code to get it. Would you please have a suggestion? Thank you so much .

paper 1 paper 2

---

***in silico* DEG details** I used intersect function in R which can do the same job. Which I am not sure is the experiment try to find transcription factors that cause the difference in phenotype between 2 conditions. DEG in bulk-RNA from DeSeq2 (log2 fold change \>2.5 or log2 fold change \< -2.5) intersect with genes differently in chromatin accessibility from DiffBind (log2 fold change \>2.5 or log2 fold change \< -2.5). I narrowed down to 142 genes.

For atac-seq result, after applying filter log2fold \> 2.5 or log2fold \< -2.5 and p value \< 0.00001, I still got around 12k genes that are different between 2 conditions. Meanwhile bulk-RNA after apply filter has only around 300 genes. The number of genes in human and pig are similar so I am not sure about the number of gene I have in atac-seq

What do you think about this approach because I don't use MEME or TOMTOM?

---

I share with you what the author did.

> " In order to screen new transcription factors by combining with star transcription factors. First, we used MEME software to predict the motif of Duroc pigs' peaks. Second, we compared the transcription factor database with TOMTOM to predict all transcription factors that may bind to the chromatin open regions . Third, we found all possible binding sites of MYOD, MYOG, MYF5, MYF6 and MEF2 families and did a union. A total of 99 sites were found. "Fourthly, we counted other transcription factors that may bind to these 99 sites and sorted them by the number of binding sites. All of the above statistics were manually counted."

Will read more to give you a good reply. Thank you so much

# Answer
You requested I make a response. Okay, I'm reluctant to do this firstly, all my eukaryotes genomics is out of date except for immunology stuff and secondly, this is a hardcore eukaryote gene expression question. I am really not sure about the pig-human side of this experiment.

What you have is an extremely good data set. What you appear to be doing is assessing the complexity of chromatin based gene expression, using very hardcore data. What we don't know is the variation e.g. in experiments leading to the RNA-seq. I would guess the RNA-seq is varying and being compared against atac-seq.

I have not read the papers, read the abstracts you posted. MEME and TOMTOM don't seem to be impressive in comparison to an absolute biological experiment RNA-seq vs. atac-seq. Predictions are never perfect (the very, very best predictions are 80%) and you have an absolute data set - why even be concerned about predictions based analysis?

I do know DEG and the approach you are using is very conservative. It's good to visualise this on a volcano plot. Personally I just use a probabilistic threshold and then assess the log2 for each gene group. You are using both -/+2.5 log2 and \<1E-6 . Personally I would be concerned with the stringency and also look at results with e.g. just \<1E-6 or \<1E-6 and -/+ 0.5log2. Keep in mind I don't do eukaryotes.

The threshold criteria will be the number 1 issue that affects your results and that is the most important information I can give. So do your approach and also try a more relaxed threshold. Personally I wouldn't relax the p\<1E-6, but I would relax the log2. The reasons for this are a bit complicated, but they are statistically sound.

Your intersect approach seems fine, viz.

```
intersect(object1, object2)

```

@haci's method is cool, the above does the job too. The question is the genes vs thresholds and the experimental set ups and then the comparison between experimental setups, i.e. which genes are switching under chromatin control.

*What I am not clear about is the experimental design and how you are performing DSeq2,*

1. I assume you are using an internal control for RNAseq and a comparable one for atac-seq. I assume you are not comparing atac-seq verse RNAseq in the DEG.
2. The alternative is atac-seq vs. RNAseq DEG and the comparison is pig versus human. I honestly dunno.

Anway the caveat I would add is an `intersect` SHOULD NOT ignore whether it is +log2 or -log2 ... you need to do four intersects

* +ve verse +ve
* -v verse -ve
* +ve verse -ve
* -ve verse +ve

Again its four separate intersects.

Good luck looks a really interesting experiment and analytics look simple but powerful in context.
> 1 votes

# Answer
If I am getting your question right, all you need is https://bioinfogp.cnb.csic.es/tools/venny/. I assume that your two lists corresponding to RNA-seq and ATAC-seq use the same set of gene identifiers (gene symbol, Entrez ID, Ensembl ID etc).

The `Results` section is populated with entries that belongs to the corresponding "region" upon clicking on the regions of the Venn diagram.

Personally I have never used it but my wet-lab colleagues use it all the time.
> 2 votes

---
Tags: rna-seq, deseq2, atac-seq
---


# Question
Title: How important are the homozygous variants that get unnecessarily deleted using liftover?
I'm referring to the text described here:

> These tools \[NCBI remap, CrossMap\] operate only on the sites present in an input VCF, and return the representation of those sites in a new genome assembly. This does not capture all variation, however. Consider an individual who has sequence reads that indicate they match the GRCh37 reference genome assembly at position GRCh37.chr1:169,519,049 (i.e. the individual's genotype is T/T). Because the individual is homozygous reference at that site, there will be no variation present in their VCF file created on GRCh37. However, the analogous position on the updated GRCh38 reference genome assembly, position GRCh38.chr1:169,549,811, has the reference base C. Consequently, if the individual's read data were analyzed on GRCh38, they would be identified as homozygous for a C-\>T SNP. Because this site is not present in the input GRCh37 VCF, it is never added when creating a GRCh38 VCF by these other tools.

How important is this? Do you think it is relevant to any types of analysis?

# Answer
> 2 votes
All this is saying is that the tool doesn't lift over variants which aren't present in the vcf. I can't imagine this being a big deal in any analysis I've done and it's exactly what I'd expect the tool to do. There's no rule that invariant sites (with respect to the reference) can't be in a vcf, so if you are concerned then just retain the variants in the vcf file at the calling stage

---
Tags: genome, assembly, liftover, ncbi, hg19
---


# Question
Title: Is a human genome file with reference as NCBI37 the same as hg19?
*This question was also asked on Biostars*

In other words, are NCBI37 and hg19 synonymous? I can't find much information for what NCBI37 is in a human (not mouse) genome.

The VCF is from a lab sequencing human genomes, and they're no longer able to be reached out to. I can't provide a link to NCBI37 because I don't know what it is and I can't find much information on it.

Relevant line from VCF: "##reference=NCBI37"

Here is part of the VCF header:

```
##fileformat=VCFv4.2
##source=apt-format-result:2.10.2
##reference=NCBI37
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=1>
##contig=<ID=10>
##contig=<ID=11>
##contig=<ID=12>
##contig=<ID=19> [omitted some chromosomes for brevity]
##contig=<ID=2>
##contig=<ID=20>
##contig=<ID=MT>
##contig=<ID=X>
##contig=<ID=Y>

```

# Answer
Most likely yes, they are the same.

However, "NCBI37" is an odd/unusual notation, because it doesn't indicate the species (human / mouse). Assuming it's a human genome (based on you comparing it to *hg*19 rather than *mm*19, or something else), then the correct name for the NCBI reference should be GRCh37.

As it's a different name, we can only make guesses as to what the people who generated the VCF meant. You'll need to talk directly to them in order to be sure.
> 1 votes

---
Tags: vcf, genome, reference-genome, hg19
---


title: Make a new column in R or bash

text:
# Question
Title: Make a new column in R or bash
I have antibody data the first column is heavy chain sequence and the second column is the light chain

I want to create a third column by combining two columns like

If `CARGRGIAVAGTTGFDYW` is a heavy chain and `CQQYNSYSLTF` a light chain, in third column I want `CARGRGIAVAGTTGFDYW_CQQYNSYSLTF`

if a heavy or light chain is absent, I want something like `NA_CQQYNSYSLTF` or `CARGRGIAVAGTTGFDYW_NA` depends of absent part

```
CDR3.light.aa   CDR3.heavy.aa
-              CAGDLSVGGVLAGGVPHLRHFDPW
-              CPETPRYCSGGFCYGEFDNW

CSSYAGSNTFGLF   -
CSSYAGSNTGLF    -

```

So far I have done this

```
a$a <- paste(a$CDR3.heavy.aa, a$CDR3.light.aa, sep= "_")

```

But I noticed a lot of mess and manual works

Any help please?

# Answer
> 2 votes
Your approach was correct, you just have to transform your dataframe first to replace the character "-" by NA :

```
library(stringr)
library(dplyr)
df=data.frame(CDR3.light.aa=c("-","-","CSSYAGSNTFGLF",'CSSYAGSNTGLF'),
              CDR3.heavy.aa=c("CAGDLSVGGVLAGGVPHLRHFDPW","-","CPETPRYCSGGFCYGEFDNW","-"))

df=df %>%
  mutate_all(funs(str_replace(., "-", "NA")))

df$a <- paste(df$CDR3.heavy.aa, df$CDR3.light.aa, sep= "_")

 df
  CDR3.light.aa            CDR3.heavy.aa                                  a
1            NA CAGDLSVGGVLAGGVPHLRHFDPW        CAGDLSVGGVLAGGVPHLRHFDPW_NA
2            NA                       NA                              NA_NA
3 CSSYAGSNTFGLF     CPETPRYCSGGFCYGEFDNW CPETPRYCSGGFCYGEFDNW_CSSYAGSNTFGLF
4  CSSYAGSNTGLF                       NA                    NA_CSSYAGSNTGLF

```

---
Tags: r
---


title: Singularity And MacOS

text:
# Question
Title: Singularity And MacOS
I am trying to install Singularity on my MacOS but it doesn't seem to pick up `which singularity` or `singularity --version`:

**Main Attempt:**

1. Install using Homebrew (and install Home-brew):

```
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"
brew install --cask --no-quarantine singularity
singularity --version

```

resulting in

```
==> Downloading https://bitbucket.org/router_gray/singularityviewer/downloads/Singularity_Alpha_1_8_9_8709_x86_64.dmg
Already downloaded: /Users/Username/Library/Caches/Homebrew/downloads/d3015dad9cc2d3ef25603a68114f3caf63cb0dd90e03f159405999262392baec--Singularity_Alpha_1_8_9_8709_x86_64.dmg
==> Uninstalling Cask singularity
==> Backing App 'SingularityAlpha.app' up to '/opt/homebrew/Caskroom/singularity/1.8.9.8709/SingularityAlpha.app'
==> Removing App '/Applications/SingularityAlpha.app'
==> Purging files for version 1.8.9.8709 of Cask singularity
==> Installing Cask singularity
==> Moving App 'SingularityAlpha.app' to '/Applications/SingularityAlpha.app'

```

Then test using `singularity --version` with the result

```
singularity not found

```

**Attempt 2:**

* Maybe it wants singularity spelled with a capital Singularity instead of lowercase:

```
which Singularity

```

resulting in:

```
Singularity not found

```

**Attempt 3:**

* Export the path to the bash\_profile and the source it

```
export PATH="/Applications/SingularityAlpha.app/Contents/MacOS:$PATH"

```

However, this resulted in the same issue when `singularity --version` was run. It could not recognize that the singularity was installed.

**Attempt 4:**

* Explore inside the Singularity app which is installed under '/Applications/SingularityAlpha.app'
* The path `/Applications/SingularityAlpha.app/Contents/MacOS` lead to a file called SingularityAlpha.
* Exporting this path in the bash profile and sourcing still lead to the same issue when trying to see if it even recoginzed the singularity using `singularity --version`.

I have also tried reinstalling it and it still has the same issue. Is there any suggestion for how to install/use Singularity locally as many packages will have to be tested locally on my machine before running it on my HPC?

# Answer
***Answer:*** Currently there is not a good way for Singularity to be installed on MacOS. Perhaps, there will be in the future.

---

I've tried and failed to install Singularity on macOS - it needs VirtualBox or some such virtualization software and a command line utility that makes virtual containers - ultimately there were so many problems that I just gave up. Docker is far easier to install on macOS and there is not a lot of difference IMO. You can use docker to build and test images and then pull them to a place where you have access to singularity and test singularity specific stuff there.

I think Singularity is a lot easier to install on Linux. See here: https://docs.sylabs.io/guides/3.0/user-guide/installation.html#mac

---

**EDIT-2**

Singularity-Alpha is a Second Life video game related cask (https://github.com/singularity-viewer/SingularityViewer), it is not the Singularity software. Don't waste your time on it, I just tried to test it and realized I'd already tested and discarded it as irrelevant.

---

EDIT: To address your homebrew issue, note that you're installing a cask and it's moving an App to /Applications/. Just typing in commands that should work is not the same as proper testing. Try opening the SingularityAlpha application to see what that does. I'm not saying it'll work, just that you can test it better. Your attempts only show a lack of understanding in how macOS Applications work and how you can(not) access them from the command line or add locations to $PATH like you can on Linux.
> 3 votes

---
Tags: singularity
---